Commit ac3dec79 authored by ymeur's avatar ymeur
Browse files

Source update to fix numerous problem for running LAM model with LMDZ physics

YM
parent 34cf5823
......@@ -3,6 +3,7 @@ MODULE time_field_mod
USE field_mod, ONLY : t_field
IMPLICIT NONE
PRIVATE
LOGICAL,PARAMETER :: xios_remap=.TRUE.
SAVE
TYPE t_time_interpolation
......@@ -17,14 +18,22 @@ MODULE time_field_mod
END type t_time_interpolation
TYPE t_time_field
CHARACTER(30) :: name, pname ! XIOS names of field and of 3D pressure used for interpolation
CHARACTER(30) :: name, pname, suffix ! XIOS names of field and of 3D pressure used for interpolation
TYPE(t_field), POINTER :: f_t(:), & ! value of field at time t in [t0,t1]
& f_t0(:), & ! value at time t0
& f_t1(:) ! value at time t1
& f_t1(:), & ! value at time t1
& f_tmp(:), &! temporary field for reading tracers
& f_levels(:) ! temporary field on read level
LOGICAL :: is_ps ! surface pressure ps plays a special role when reading 3D input files
LOGICAL :: is_q !
END type t_time_field
TYPE(t_field), POINTER :: f_p(:) ! 3D pressure field reconstructed from ps and hybrid coordinate
TYPE(t_field), POINTER :: f_p(:) ! 3D pressure field reconstructed from ps and hybrid coordinate
TYPE(t_field), POINTER :: f_ps(:) ! ps
REAL(rstd),ALLOCATABLE :: levels(:)
INTEGER :: nlevels
PUBLIC :: t_time_interpolation, init_time_interpolation, update_time_interpolation, &
& t_time_field, init_time_field, allocate_time_field, &
......@@ -33,9 +42,23 @@ MODULE time_field_mod
CONTAINS
SUBROUTINE init_time_field
USE icosa
USE field_mod, ONLY : field_T, type_real, allocate_field
USE grid_param, ONLY : llm
USE xios_mod
USE omp_para
CALL allocate_field(f_p, field_t, type_real, llm, name='pressure_time_field')
CALL allocate_field(f_ps, field_t, type_real, name='pressure_surface_time_field')
IF (is_omp_master) CALL xios_get_axis_attr("axis_nudged",n_glo=nlevels)
CALL bcast_omp(nlevels)
ALLOCATE(levels(nlevels))
IF (is_omp_master) CALL xios_get_axis_attr("axis_nudged",value=levels)
CALL bcast_omp(levels)
! levels=levels*100 ! millibars -> Pa
END SUBROUTINE init_time_field
SUBROUTINE init_time_interpolation(interval, interp)
......@@ -44,30 +67,41 @@ CONTAINS
interp%interval = interval
END SUBROUTINE init_time_interpolation
SUBROUTINE allocate_time_field(field, field_type, data_type, dim3, dim4, name, pname)
SUBROUTINE allocate_time_field(field, field_type, data_type, dim3, dim4, name,suffix, pname)
USE icosa
USE field_mod, ONLY : allocate_field
USE domain_mod, ONLY : ndomain, assigned_domain
USE xios_mod
USE omp_para
! value of argument pname must reflect the contents of XIOS XML configuration files
TYPE(t_time_field), INTENT(OUT) :: field
INTEGER, INTENT(IN) :: field_type
INTEGER, INTENT(IN) :: data_type
INTEGER, OPTIONAL, INTENT(IN) :: dim3, dim4
CHARACTER(*), OPTIONAL, INTENT(IN) :: name, pname
CHARACTER(*), OPTIONAL, INTENT(IN) :: name, suffix, pname
INTEGER :: ind
field%name = TRIM(name) ! note that name is in fact a compulsory argument
field%suffix =""
IF (PRESENT(suffix)) field%suffix = TRIM(suffix)
field%is_ps = PRESENT(pname)
IF(field%is_ps) field%pname = TRIM(pname)
field%is_q = .FALSE.
IF (PRESENT(dim4) .AND. dim4 == nqtot) field%is_q = .TRUE.
CALL allocate_field(field%f_t, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t')
CALL allocate_field(field%f_t0, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t0')
CALL allocate_field(field%f_t1, field_type, data_type, dim3, dim4, name=TRIM(name)//'_t1')
CALL allocate_field(field%f_tmp, field_type, data_type, dim3, name=TRIM(name)//'_tmp')
CALL allocate_field(field%f_levels, field_type, data_type, nlevels, name=TRIM(name)//'_level')
DO ind=1, ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
field%f_t(ind)%rval4d = 0.
field%f_t0(ind)%rval4d = 0.
field%f_t1(ind)%rval4d = 0.
field%f_tmp(ind)%rval4d = 0.
field%f_levels(ind)%rval4d = 0.
END DO
!$OMP BARRIER
END SUBROUTINE allocate_time_field
......@@ -135,15 +169,46 @@ CONTAINS
USE layout_mod, ONLY : field_T, field_U, field_Z
USE transfert_mod, ONLY : transfert_request, req_i0, req_i1
USE pression_mod, ONLY : pression_mid
USE xios_mod, ONLY : xios_read_field, xios_write_field
USE xios_mod, ONLY : xios_read_field, xios_write_field, xios_is_valid_field
USE tracer_icosa_mod, ONLY : tracers
USE icosa
USE write_field_mod
USE omp_para
USE vertical_remap_mod
TYPE(t_time_interpolation), INTENT(IN) :: interp
TYPE(t_time_field), INTENT(INOUT) :: field
REAL(rstd) :: t0,t1
INTEGER :: iq,l,ind
IF(interp%read_t0) THEN
! we get data for t0 only once, from XIOS field XXX0
CALL xios_read_field(TRIM(field%name)//"0", field%f_t0)
IF (field%is_q) THEN ! for tracers
DO iq=1,nqtot
! CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix)//"0", field%f_tmp)
IF (xios_is_valid_field(TRIM(tracers(iq)%name)//TRIM(field%suffix))) THEN
CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix)//"0", field%f_tmp)
ELSE
DO ind=1,ndomain
field%f_tmp(ind)%rval4d(:,:,:)=0
ENDDO
ENDIF
DO ind=1,ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
CALL swap_dimensions(ind)
CALL swap_geometry(ind)
DO l=ll_begin,ll_end
field%f_t0(ind)%rval4d(:,l,iq) = field%f_tmp(ind)%rval4d(:,l,1)
ENDDO
ENDDO
ENDDO
ELSE
CALL xios_read_field(TRIM(field%name)//"0", field%f_t0)
! CALL writeField(TRIM(field%name),field%f_t0)
ENDIF
SELECT CASE(field%f_t0(1)%field_type)
CASE(field_T)
CALL transfert_request(field%f_t0, req_i0)
......@@ -153,15 +218,44 @@ CONTAINS
END SELECT
IF(field%is_ps) THEN
DO ind=1,ndomain
f_ps(ind)%rval4d = field%f_t0(ind)%rval4d
ENDDO
CALL pression_mid(field%f_t0, f_p)
CALL xios_write_field(TRIM(field%pname)//"0",f_p)
END IF
! CALL writeField(TRIM(field%pname),f_p)
END IF
END IF
IF(interp%swap) CALL swap_time_field(field)
IF(interp%read_t1) THEN
CALL xios_read_field(TRIM(field%name), field%f_t1)
IF (field%is_q) THEN ! for tracers
DO iq=1,nqtot
IF (xios_is_valid_field(TRIM(tracers(iq)%name)//TRIM(field%suffix))) THEN
CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix), field%f_tmp)
ELSE
DO ind=1,ndomain
field%f_tmp(ind)%rval4d(:,:,:)=0
ENDDO
ENDIF
DO ind=1,ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
CALL swap_dimensions(ind)
CALL swap_geometry(ind)
DO l=ll_begin,ll_end
field%f_t1(ind)%rval4d(:,l,iq) = field%f_tmp(ind)%rval4d(:,l,1)
ENDDO
ENDDO
ENDDO
ELSE
CALL xios_read_field(TRIM(field%name), field%f_t1)
! CALL writeField(TRIM(field%name),field%f_t1)
ENDIF
SELECT CASE(field%f_t1(1)%field_type)
CASE(field_T)
CALL transfert_request(field%f_t1, req_i0)
......@@ -173,12 +267,144 @@ CONTAINS
IF(field%is_ps) THEN
CALL pression_mid(field%f_t1, f_p)
CALL xios_write_field(field%pname, f_p)
! CALL writeField(TRIM(field%pname),f_p)
END IF
ENDIF
END SUBROUTINE read_time_field
SUBROUTINE read_time_field_remap(interp, field)
! surface pressure ps must always be read before 3D fields
! in the special case that we read ps, we compute the 3D pressure field
! and provide it to XIOS for later use in vertical interpolation of 3D fields
! this special role of ps is marked by the flag is_ps
USE layout_mod, ONLY : field_T, field_U, field_Z
USE transfert_mod, ONLY : transfert_request, req_i0, req_i1
USE pression_mod, ONLY : pression_mid
USE xios_mod, ONLY : xios_read_field, xios_write_field, xios_is_valid_field
USE tracer_icosa_mod, ONLY : tracers
USE icosa
USE write_field_mod
USE omp_para
USE vertical_remap_mod
TYPE(t_time_interpolation), INTENT(IN) :: interp
TYPE(t_time_field), INTENT(INOUT) :: field
REAL(rstd) :: t0,t1
INTEGER :: iq,l,ind
IF(interp%read_t0) THEN
! we get data for t0 only once, from XIOS field XXX0
IF (field%is_q) THEN ! for tracers
DO iq=1,nqtot
! CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix)//"0", field%f_tmp)
IF (xios_is_valid_field(TRIM(tracers(iq)%name)//TRIM(field%suffix))) THEN
CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix)//"0", field%f_levels)
ELSE
DO ind=1,ndomain
field%f_levels(ind)%rval4d(:,:,:)=0
ENDDO
ENDIF
CALL vertical_remap(levels,field%f_levels,f_ps,field%f_tmp)
DO ind=1,ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
CALL swap_dimensions(ind)
CALL swap_geometry(ind)
DO l=ll_begin,ll_end
field%f_t0(ind)%rval4d(:,l,iq) = field%f_tmp(ind)%rval4d(:,l,1)
ENDDO
ENDDO
ENDDO
ELSE
IF(field%is_ps) THEN
CALL xios_read_field(TRIM(field%name)//"0", field%f_t0)
ELSE
CALL xios_read_field(TRIM(field%name)//"0", field%f_levels)
CALL vertical_remap(levels,field%f_levels,f_ps,field%f_t0)
ENDIF
CALL writeField(TRIM(field%name),field%f_t0)
ENDIF
SELECT CASE(field%f_t0(1)%field_type)
CASE(field_T)
CALL transfert_request(field%f_t0, req_i0)
CALL transfert_request(field%f_t0, req_i1)
CASE DEFAULT
STOP 'read_time_field : unsupported field_type'
END SELECT
IF(field%is_ps) THEN
DO ind=1,ndomain
f_ps(ind)%rval4d = field%f_t0(ind)%rval4d
ENDDO
CALL pression_mid(field%f_t0, f_p)
CALL xios_write_field(TRIM(field%pname)//"0",f_p)
CALL writeField(TRIM(field%pname),f_p)
END IF
END IF
IF(interp%swap) CALL swap_time_field(field)
IF(interp%read_t1) THEN
IF (field%is_q) THEN ! for tracers
DO iq=1,nqtot
! CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix), field%f_tmp)
IF (xios_is_valid_field(TRIM(tracers(iq)%name)//TRIM(field%suffix))) THEN
CALL xios_read_field(TRIM(tracers(iq)%name)//TRIM(field%suffix), field%f_levels)
ELSE
DO ind=1,ndomain
field%f_levels(ind)%rval4d(:,:,:)=0
ENDDO
ENDIF
CALL vertical_remap(levels,field%f_levels,f_ps,field%f_tmp)
DO ind=1,ndomain
IF(.NOT. assigned_domain(ind)) CYCLE
CALL swap_dimensions(ind)
CALL swap_geometry(ind)
DO l=ll_begin,ll_end
field%f_t1(ind)%rval4d(:,l,iq) = field%f_tmp(ind)%rval4d(:,l,1)
ENDDO
ENDDO
ENDDO
ELSE
IF(field%is_ps) THEN
CALL xios_read_field(TRIM(field%name), field%f_t1)
ELSE
CALL xios_read_field(TRIM(field%name), field%f_levels)
CALL vertical_remap(levels,field%f_levels,f_ps,field%f_t1)
ENDIF
! CALL writeField(TRIM(field%name),field%f_t1)
ENDIF
SELECT CASE(field%f_t1(1)%field_type)
CASE(field_T)
CALL transfert_request(field%f_t1, req_i0)
CALL transfert_request(field%f_t1, req_i1)
CASE DEFAULT
STOP 'read_time_field : unsupported field_type'
END SELECT
IF(field%is_ps) THEN
CALL pression_mid(field%f_t1, f_p)
CALL xios_write_field(field%pname, f_p)
! CALL writeField(TRIM(field%pname),f_p)
END IF
ENDIF
END SUBROUTINE read_time_field_remap
SUBROUTINE interpolate_time_field(interp, field)
USE field_mod, ONLY : t_field
USE domain_mod, ONLY : ndomain, assigned_domain
......
......@@ -69,8 +69,8 @@ CONTAINS
TYPE(t_field),POINTER :: f_theta_rhodz(:)
TYPE(t_field),POINTER :: f_u(:)
TYPE(t_field),POINTER :: f_q(:)
IF (metric_type == metric_icosa_area) CALL up_border(tt, f_ps, f_mass, f_theta_rhodz, f_u, f_q)
IF (metric_type == metric_icosa_area ) CALL up_border(tt, f_ps, f_mass, f_theta_rhodz, f_u, f_q)
END SUBROUTINE update_border
......
......@@ -5,7 +5,7 @@ MODULE nudging_mod
PRIVATE
SAVE
INTEGER,PARAMETER :: halo_scheme=4
INTEGER,PARAMETER,PUBLIC :: lam_halo_scheme=5
TYPE(t_time_interpolation) :: interp
!$OMP THREADPRIVATE(interp)
......@@ -94,7 +94,7 @@ CONTAINS
CALL allocate_time_field(tf_ulat, field_T, type_real, llm, name='ulat_nudged')
CALL allocate_time_field(tf_u, field_U, type_real, llm, name='u_nudged')
CALL allocate_time_field(tf_temp, field_t, type_real, llm, name='temp_nudged')
CALL allocate_time_field(tf_q, field_T, type_real, llm, nqtot, name='q_nudged')
CALL allocate_time_field(tf_q, field_T, type_real, llm, nqtot, name='q', suffix="_nudged")
CALL allocate_time_field(tf_theta_rhodz, field_T, type_real, llm, 1, name='theta_rhodz_nudged')
CALL getin("guide_T",guide_T)
......@@ -551,7 +551,7 @@ CONTAINS
ENDIF
!$OMP BARRIER
DO n=1, NINT(halo_scheme+stiffness)
DO n=1, NINT(lam_halo_scheme+stiffness)
IF (is_omp_master) THEN
DO ind = 1 , ndomain
CALL swap_dimensions(ind)
......@@ -594,11 +594,11 @@ CONTAINS
IF (neighbours(ij)==-1) THEN
coeff_i(ij)=0
border(ij) = 1
ELSE IF (neighbours(ij)<=halo_scheme) THEN
ELSE IF (neighbours(ij)<=lam_halo_scheme) THEN
coeff_i(ij)=1
border(ij) = 0
ELSE
coeff_i(ij) = 0.5*tanh(-2.5*(neighbours(ij)-halo_scheme-stiffness/2)/stiffness)+0.5
coeff_i(ij) = 0.5*tanh(-2.5*(neighbours(ij)-lam_halo_scheme-stiffness/2)/stiffness)+0.5
border(ij) = 1
ENDIF
ENDDO
......@@ -839,6 +839,7 @@ CONTAINS
SUBROUTINE guided(time, f_ps, f_mass, f_theta_rhodz, f_u, f_q)
USE disvert_mod
USE checksum_mod
USE write_field_mod
IMPLICIT NONE
REAL(rstd), INTENT(IN) :: time
TYPE(t_field), POINTER :: f_u(:)! initial condition
......@@ -861,7 +862,7 @@ CONTAINS
INTEGER :: ind
CALL check_time_interp(time)
!$OMP BARRIER
DO ind = 1 , ndomain
IF (.NOT. assigned_domain(ind)) CYCLE
......
......@@ -213,7 +213,7 @@ CONTAINS
CALL MPI_ALLGATHERV(ind_glo, ncell, MPI_INTEGER, ind_glo_tot, ncell_mpi, displ, MPI_INTEGER, comm_icosa,ierr)
DO i=1,ncell_tot
cell_glo_tot(ind_glo_tot(i))= 0
cell_glo_tot(ind_glo_tot(i))= cell_glo_tot(ind_glo_tot(i)) + 1
ENDDO
ncell_glo=0
......@@ -229,7 +229,7 @@ CONTAINS
ENDDO
CALL xios_set_domain_attr("u",ni_glo=ncell_tot, ibegin=displ(mpi_rank), ni=ncell)
CALL xios_set_domain_attr("u",ni_glo=ncell_glo, ibegin=displ(mpi_rank), ni=ncell)
CALL xios_set_domain_attr("u", data_dim=1, type='unstructured' , nvertex=2, i_index=ind_glo)
CALL xios_set_domain_attr("u",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
......
......@@ -293,6 +293,7 @@ CONTAINS
ENDDO
! assign edge and vertex in the inner domain
nbp=0
DO ind_d=1,ndomain_glo
d=>domain_glo(ind_d)
......@@ -389,6 +390,79 @@ CONTAINS
! assign missing cells and edge, vertex for limited area config
! first assign cells and vertex on the first neighbourough
DO ind_d=1,ndomain_glo
d=>domain_glo(ind_d)
nf=d%face
DO j=d%jj_begin,d%jj_end
DO i=d%ii_begin,d%ii_end
ii=d%ii_begin_glo-d%ii_begin+i
jj=d%jj_begin_glo-d%jj_begin+j
ind=vertex_glo(ii,jj,nf)%ind
delta=vertex_glo(ii,jj,nf)%delta
IF (cell_glo(ind)%assign_face==nf) THEN
IF (cell_glo(ind)%assign_domain==-1) THEN
cell_glo(ind)%assign_domain=ind_d
cell_glo(ind)%assign_i=i
cell_glo(ind)%assign_j=j
ENDIF
ENDIF
CALL assign_edge_if_not(ind_d,ind,i,j,delta,0)
CALL assign_edge_if_not(ind_d,ind,i,j,delta,1)
CALL assign_edge_if_not(ind_d,ind,i,j,delta,2)
CALL assign_edge_if_not(ind_d,ind,i,j,delta,3)
CALL assign_edge_if_not(ind_d,ind,i,j,delta,4)
CALL assign_edge_if_not(ind_d,ind,i,j,delta,5)
CALL assign_vertex_if_not(ind_d,ind,i,j,delta,0)
CALL assign_vertex_if_not(ind_d,ind,i,j,delta,1)
CALL assign_vertex_if_not(ind_d,ind,i,j,delta,2)
CALL assign_vertex_if_not(ind_d,ind,i,j,delta,3)
CALL assign_vertex_if_not(ind_d,ind,i,j,delta,4)
CALL assign_vertex_if_not(ind_d,ind,i,j,delta,5)
ENDDO
ENDDO
ENDDO
! assign edges and vertex on the second neighbourough
DO ind_d=1,ndomain_glo
d=>domain_glo(ind_d)
nf=d%face
DO j=d%jj_begin-1,d%jj_end+1
DO i=d%ii_begin-1,d%ii_end+1
ii=d%ii_begin_glo-d%ii_begin+i
jj=d%jj_begin_glo-d%jj_begin+j
ind=vertex_glo(ii,jj,nf)%ind
delta=vertex_glo(ii,jj,nf)%delta
IF (cell_glo(ind)%assign_face==nf) THEN
IF (cell_glo(ind)%assign_domain==-1) THEN
cell_glo(ind)%assign_domain=ind_d
cell_glo(ind)%assign_i=i
cell_glo(ind)%assign_j=j
ENDIF
ENDIF
IF ( i+1 <= d%ii_end+1 ) CALL assign_edge_if_not(ind_d,ind,i,j,delta,0)
IF ( j+1 <= d%jj_end+1 ) CALL assign_edge_if_not(ind_d,ind,i,j,delta,1)
IF ( i-1 >= d%ii_begin-1 .AND. j+1<=d%jj_end+1 ) CALL assign_edge_if_not(ind_d,ind,i,j,delta,2)
IF ( i-1 >= d%ii_begin-1 ) CALL assign_edge_if_not(ind_d,ind,i,j,delta,3)
IF ( j-1 >= d%jj_begin-1 ) CALL assign_edge_if_not(ind_d,ind,i,j,delta,4)
IF ( i+1 <= d%ii_end+1 .AND. j-1 >=d%jj_begin-1 ) CALL assign_edge_if_not(ind_d,ind,i,j,delta,5)
IF ( i+1 <= d%ii_end+1 .AND. j+1 <= d%jj_end+1) CALL assign_vertex_if_not(ind_d,ind,i,j,delta,0)
IF ( i-1 >= d%ii_begin-1 .AND. j+1 <= d%jj_end+1) CALL assign_vertex_if_not(ind_d,ind,i,j,delta,1)
IF ( i-1 >= d%ii_begin-1 .AND. j+1 <= d%jj_end+1) CALL assign_vertex_if_not(ind_d,ind,i,j,delta,2)
IF ( i-1 >= d%ii_begin-1 .AND. j-1 >= d%jj_begin-1) CALL assign_vertex_if_not(ind_d,ind,i,j,delta,3)
IF ( i+1 <= d%ii_end+1 .AND. j-1 >= d%jj_begin-1) CALL assign_vertex_if_not(ind_d,ind,i,j,delta,4)
IF ( i+1 <= d%ii_end+1 .AND. j-1 >= d%jj_begin-1) CALL assign_vertex_if_not(ind_d,ind,i,j,delta,5)
ENDDO
ENDDO
ENDDO
! assign edges and vertex on the third neighbourough
DO ind_d=1,ndomain_glo
d=>domain_glo(ind_d)
nf=d%face
......@@ -422,7 +496,8 @@ CONTAINS
ENDDO
ENDDO
ENDDO
ENDDO
DO ind_d=1,ndomain_glo
......@@ -481,6 +556,19 @@ CONTAINS
ENDDO
ENDDO
nbp=0
DO ind_d=1,ndomain_glo
d=>domain_glo(ind_d)
nf=d%face
DO j=d%jj_begin,d%jj_end
DO i=d%ii_begin,d%ii_end
DO k=0,5
IF (d%edge_assign_domain(k,i,j)==ind_d) nbp=nbp+1
ENDDO
ENDDO
ENDDO
ENDDO
CONTAINS
SUBROUTINE assign_edge(ind_d,ind,i,j,delta,k)
......
......@@ -372,6 +372,7 @@ CONTAINS
CALL enter_profile(id_adv)
IF(MOD(it,itau_adv)==0) THEN
CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt) ! accumulate mass and fluxes if diagflux_on
fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
......@@ -396,6 +397,7 @@ CONTAINS
IF(positive_theta) CALL abort_acc("positive_theta")
IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
END IF
CALL exit_profile(id_adv)
! CALL checksum(f_ps)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment