MODULE interface_icosa_lmdz_mod USE field_mod, ONLY: t_field USE transfert_mod, ONLY: t_message TYPE(t_message),SAVE :: req_u, req_z TYPE(t_message),SAVE :: req_dps0, req_dulon0, req_dulat0, req_dTemp0, req_dq0 TYPE(t_field),POINTER,SAVE :: f_p(:) TYPE(t_field),POINTER,SAVE :: f_pks(:) TYPE(t_field),POINTER,SAVE :: f_pk(:) TYPE(t_field),POINTER,SAVE :: f_p_layer(:) TYPE(t_field),POINTER,SAVE :: f_theta(:) TYPE(t_field),POINTER,SAVE :: f_phi(:) TYPE(t_field),POINTER,SAVE :: f_Temp(:) TYPE(t_field),POINTER,SAVE :: f_ulon(:) TYPE(t_field),POINTER,SAVE :: f_ulat(:) TYPE(t_field),POINTER,SAVE :: f_vort(:) TYPE(t_field),POINTER,SAVE :: f_vortc(:) TYPE(t_field),POINTER,SAVE :: f_dulon(:) TYPE(t_field),POINTER,SAVE :: f_dulat(:) TYPE(t_field),POINTER,SAVE :: f_dTemp(:) TYPE(t_field),POINTER,SAVE :: f_dq(:) TYPE(t_field),POINTER,SAVE :: f_dps(:) TYPE(t_field),POINTER,SAVE :: f_duc(:) TYPE(t_field),POINTER,SAVE :: f_bounds_lon(:) TYPE(t_field),POINTER,SAVE :: f_bounds_lat(:) INTEGER :: start_clock INTEGER :: stop_clock INTEGER :: count_clock=0 INTEGER,SAVE :: nbp_phys INTEGER,SAVE :: nbp_phys_glo REAL,ALLOCATABLE,SAVE :: q_ave(:, :) !$OMP THREADPRIVATE(q_ave) REAL, PARAMETER :: q_epsilon = 1e-8 CONTAINS SUBROUTINE pre_initialize_physics USE etat0_plugin_mod USE icolmdz_etat0, ONLY: init_etat0_lmdz => init_etat0, etat0_lmdz => etat0 USE icolmdz_param_gravity_wave, ONLY: init_param_gravity_wave USE isotopes_mod, ONLY : using_iso IMPLICIT NONE !$OMP PARALLEL init_etat0_plugin => init_etat0_lmdz etat0_plugin => etat0_lmdz CALL init_param_gravity_wave !$OMP END PARALLEL END SUBROUTINE pre_initialize_physics SUBROUTINE initialize_physics USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz ! from dynamico USE domain_mod USE dimensions USE mpi_mod USE mpipara USE disvert_mod USE xios_mod USE time_mod , init_time_icosa=> init_time USE transfert_mod USE nudging_mod, ONLY : lam_halo_scheme ! from LMDZ USE mod_grid_phy_lmdz, ONLY : unstructured USE mod_phys_lmdz_omp_data, ONLY: klon_omp USE transfert_mod USE physics_distribution_mod, ONLY : init_physics_distribution USE infotrac_phy, ONLY: init_infotrac_phy USE icolmdz_param_gravity_wave, ONLY: param_gravity_wave IMPLICIT NONE INTEGER :: ind,i,j,ij,pos,h REAL(rstd),POINTER :: bounds_lon(:,:) REAL(rstd),POINTER :: bounds_lat(:,:) REAL(rstd),ALLOCATABLE :: latfi(:) REAL(rstd),ALLOCATABLE :: lonfi(:) REAL(rstd),ALLOCATABLE :: airefi(:) REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:) REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:) LOGICAL ,ALLOCATABLE :: outside(:,:) LOGICAL ,ALLOCATABLE :: outside_tmp(:,:) LOGICAL ,POINTER :: out(:,:) ! REAL(rstd) :: pseudoalt(llm) INTEGER :: nbp_phys, nbp_phys_glo !$OMP PARALLEL CALL allocate_field(f_bounds_lon,field_t,type_real,6) CALL allocate_field(f_bounds_lat,field_t,type_real,6) CALL allocate_field(f_p,field_t,type_real,llm+1,name="p_in") CALL allocate_field(f_pks,field_t,type_real) CALL allocate_field(f_pk,field_t,type_real,llm) CALL allocate_field(f_p_layer,field_t,type_real,llm,name="p_layer_in") CALL allocate_field(f_theta,field_t,type_real,llm) CALL allocate_field(f_phi,field_t,type_real,llm,name="phi_in") CALL allocate_field(f_Temp,field_t,type_real,llm,name="Temp_in") CALL allocate_field(f_ulon,field_t,type_real,llm,name="ulon_in") CALL allocate_field(f_ulat,field_t,type_real,llm,name="ulat_in") CALL allocate_field(f_vort,field_z,type_real,llm,name="vort_in") CALL allocate_field(f_vortc,field_t,type_real,llm,name="vortc_in") CALL allocate_field(f_dulon,field_t,type_real,llm,name="dulon_out") CALL allocate_field(f_dulat,field_t,type_real,llm,name="dulat_out") CALL allocate_field(f_dTemp,field_t,type_real,llm,name="dTemp_out") CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,name="dq_out") CALL allocate_field(f_dps,field_t,type_real,name="dps_out") CALL allocate_field(f_duc,field_t,type_real,3,llm) CALL init_message(f_dps,req_i0,req_dps0) CALL init_message(f_dulon,req_i0,req_dulon0) CALL init_message(f_dulat,req_i0,req_dulat0) CALL init_message(f_dTemp,req_i0,req_dTemp0) CALL init_message(f_dq,req_i0,req_dq0) !$OMP END PARALLEL nbp_phys=0 DO ind=1,ndomain CALL swap_dimensions(ind) ALLOCATE(outside(ii_begin:ii_end,jj_begin:jj_end)) ! for limited area : don't take cells arround the border ALLOCATE(outside_tmp(ii_begin-1:ii_end+1,jj_begin-1:jj_end+1)) ! for limited area : don't take cells arround the border out=>domain(ind)%outside DO j=jj_begin,jj_end DO i=ii_begin,ii_end outside(i,j)= out(i+1,j) .OR. & ! right out(i,j+1 ) .OR. & ! rup out(i-1 ,j+1) .OR. & ! lup out(i-1 ,j) .OR. & !left out(i ,j-1) .OR. & !ldown out(i+1,j-1) !rdown ENDDO ENDDO outside_tmp=.FALSE. outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside DO h=1,lam_halo_scheme-1 ! do not compute physic on limited area halo DO j=jj_begin,jj_end DO i=ii_begin,ii_end outside(i,j) = outside_tmp(i,j) .OR. & outside_tmp(i+1,j) .OR. & ! right outside_tmp(i,j+1 ) .OR. & ! rup outside_tmp(i-1 ,j+1) .OR. & ! lup outside_tmp(i-1 ,j) .OR. & !left outside_tmp(i ,j-1) .OR. & !ldown outside_tmp(i+1,j-1) !rdown ENDDO ENDDO outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside ENDDO DO j=jj_begin,jj_end DO i=ii_begin,ii_end IF (domain(ind)%own(i,j) .AND. .NOT.outside(i,j)) nbp_phys=nbp_phys+1 ENDDO ENDDO DEALLOCATE(outside) DEALLOCATE(outside_tmp) ENDDO !initialize LMDZ5 physic mpi decomposition CALL MPI_ALLREDUCE(nbp_phys,nbp_phys_glo,1,MPI_INTEGER,MPI_SUM,comm_icosa,ierr) CALL init_physics_distribution(unstructured, 6, nbp_phys, 1, nbp_phys_glo, llm, comm_icosa) DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) bounds_lon=f_bounds_lon(ind) bounds_lat=f_bounds_lat(ind) DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i CALL xyz2lonlat(xyz_v(ij+z_rup,:), bounds_lon(ij,1), bounds_lat(ij,1)) CALL xyz2lonlat(xyz_v(ij+z_up,:), bounds_lon(ij,2), bounds_lat(ij,2)) CALL xyz2lonlat(xyz_v(ij+z_lup,:), bounds_lon(ij,3), bounds_lat(ij,3)) CALL xyz2lonlat(xyz_v(ij+z_ldown,:), bounds_lon(ij,4), bounds_lat(ij,4)) CALL xyz2lonlat(xyz_v(ij+z_down,:), bounds_lon(ij,5), bounds_lat(ij,5)) CALL xyz2lonlat(xyz_v(ij+z_rdown,:), bounds_lon(ij,6), bounds_lat(ij,6)) ENDDO ENDDO ENDDO CALL init_infotrac_phy !$OMP PARALLEL CALL initialize_physics_omp CALL param_gravity_wave !$OMP END PARALLEL CALL xios_set_context END SUBROUTINE initialize_physics SUBROUTINE initialize_physics_omp USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz ! from dynamico USE domain_mod USE dimensions USE mpi_mod USE mpipara USE disvert_mod USE earth_const, ONLY: scale_height USE xios_mod USE time_mod , init_time_icosa=> init_time USE omp_para ! from LMDZ USE mod_grid_phy_lmdz, ONLY : unstructured, klon_glo USE mod_phys_lmdz_para, ONLY: klon_omp, reduce_min_lmdz => reduce_min , gather_lmdz => gather , bcast_lmdz => bcast USE time_phylmdz_mod, ONLY: init_time_lmdz => init_time USE transfert_mod ! USE physics_distribution_mod, ONLY : init_physics_distribution USE geometry_mod, ONLY : init_geometry USE vertical_layers_mod, ONLY : init_vertical_layers USE infotrac_phy, ONLY : init_infotrac_phy USE inifis_mod, ONLY : inifis USE readTracFiles_mod, ONLY: trac_type, isot_type USE tracer_icosa_mod, ONLY : tracs USE readTracFiles_mod, ONLY: delPhase ! USE phyaqua_mod, ONLY : iniaqua IMPLICIT NONE INTEGER :: ind,i,j,k,ij,pos REAL(rstd),POINTER :: bounds_lon(:,:) REAL(rstd),POINTER :: bounds_lat(:,:) REAL(rstd),ALLOCATABLE :: latfi(:) REAL(rstd),ALLOCATABLE :: lonfi(:) REAL(rstd),ALLOCATABLE :: airefi(:) REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:) REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:) REAL(rstd),ALLOCATABLE :: ind_cell_glo_r(:) INTEGER, ALLOCATABLE :: ind_cell_glo(:) INTEGER, ALLOCATABLE :: ind_cell_glo_tot(:) INTEGER, ALLOCATABLE :: cell_glo_tot(:) INTEGER :: ncell_glo_tot REAL(rstd) :: pseudoalt(llm) REAL(rstd) :: aps(llm) REAL(rstd) :: bps(llm) REAL(rstd) :: scaleheight INTEGER :: run_length REAL :: day_length ! length of a day (s) ! SAVEd to be OpenMP shared <--- NO!!!! INTEGER :: annee_ref INTEGER :: day_ref INTEGER :: day_ini REAL :: start_time REAL :: physics_timestep ! Tracer stuff (SAVEd when needed to be OpenMP shared) INTEGER :: nq INTEGER :: nqo, nbtr, nbtr_inca CHARACTER(len=256) :: type_trac INTEGER,ALLOCATABLE :: conv_flg(:) ! conv_flg(it)=0 : convection desactivated for tracer number it INTEGER,ALLOCATABLE :: pbl_flg(:) ! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it CHARACTER(len=8),ALLOCATABLE :: solsym(:) ! tracer name from inca TYPE(t_field),POINTER,SAVE :: f_ind_cell_glo(:) INTEGER :: iflag_phys INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca ! index of horizontal trasport schema INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca ! index of vertical trasport schema TYPE(trac_type) :: tracers_ico2lmdz(nqtot) !=== TRACERS DESCRIPTORS VECTOR TYPE(isot_type) :: isotopes_ico2lmdz(1) !=== ISOTOPES PARAMETERS VECTOR INTEGER :: iq CHARACTER(LEN=3) :: descrq(30) !--- Advection scheme description tags logical, save :: first = .TRUE. CALL init_distrib_icosa_lmdz ALLOCATE(latfi(klon_omp)) ALLOCATE(lonfi(klon_omp)) ALLOCATE(airefi(klon_omp)) ALLOCATE(bounds_latfi(klon_omp,6)) ALLOCATE(bounds_lonfi(klon_omp,6)) ALLOCATE(ind_cell_glo_r(klon_omp)) ALLOCATE(ind_cell_glo(klon_omp)) ALLOCATE(q_ave(llm,nqtot)) CALL transfer_icosa_to_lmdz(geom%lat_i,latfi) CALL transfer_icosa_to_lmdz(geom%lon_i,lonfi) CALL transfer_icosa_to_lmdz(f_bounds_lat,bounds_latfi) CALL transfer_icosa_to_lmdz(f_bounds_lon,bounds_lonfi) CALL transfer_icosa_to_lmdz(geom%Ai,airefi) CALL allocate_field(f_ind_cell_glo,field_t,type_real) DO ind=1,ndomain IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE CALL swap_dimensions(ind) CALL swap_geometry(ind) DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i f_ind_cell_glo(ind)%rval2d(ij)=domain(ind)%assign_cell_glo(i,j) ENDDO ENDDO ENDDO CALL transfer_icosa_to_lmdz(f_ind_cell_glo,ind_cell_glo_r) CALL deallocate_field(f_ind_cell_glo) ind_cell_glo=INT(ind_cell_glo_r) DEALLOCATE(ind_cell_glo_r) CALL reduce_min_lmdz(MINVAL(-ind_cell_glo),ncell_glo_tot) ! reduce_max does not exist in lmdz, use reduce_min CALL bcast_lmdz(ncell_glo_tot) ncell_glo_tot=-ncell_glo_tot ALLOCATE(cell_glo_tot(0:ncell_glo_tot)) ALLOCATE(ind_cell_glo_tot(klon_glo)) CALL gather_lmdz(ind_cell_glo,ind_cell_glo_tot) CALL bcast_lmdz(ind_cell_glo_tot) cell_glo_tot=-1 DO i=1,klon_glo cell_glo_tot(ind_cell_glo_tot(i))= 0 ENDDO pos=0 DO i=0,ncell_glo_tot IF (cell_glo_tot(i)/=-1) THEN cell_glo_tot(i) = pos pos=pos + 1 ENDIF ENDDO DO i=1,klon_omp ind_cell_glo(i)=cell_glo_tot(ind_cell_glo(i)) ENDDO ind_cell_glo = ind_cell_glo + 1 ! lmdz expect global indices begining to 1 not 0 ! CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_mpi,1,MPI_INTEGER,comm_icosa,ierr) ! ! displ(0)=0 ! DO i=1,mpi_size-1 ! displ(i)=displ(i-1)+ncell_mpi(i-1) ! ENDDO ! ! ALLOCATE(ind_glo_tot(ncell_tot)) ! ALLOCATE(cell_glo_tot(0:ncell_glo_tot-1)) ! ! cell_glo_tot(:)= -1 ! 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 ! ENDDO ! ! ncell_glo=0 ! DO i=0,ncell_glo_tot-1 ! IF (cell_glo_tot(i)/=-1) THEN ! cell_glo_tot(i) = ncell_glo ! ncell_glo=ncell_glo + 1 ! ENDIF ! ENDDO ! ! DO i=1,ncell ! ind_glo(i)=cell_glo_tot(ind_glo(i)) ! ENDDO CALL init_geometry(klon_omp,lonfi, latfi, bounds_lonfi, bounds_latfi, airefi, ind_cell_glo) scaleheight=scale_height/1000. ! Atmospheric scale height (km) aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1)) bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1)) pseudoalt(:)=-scaleheight*log(presnivs(:)/preff) CALL init_vertical_layers(llm,preff,scaleheight,ap,bp,aps,bps,presnivs,presinter,pseudoalt) ! Initialize tracer names, numbers, etc. for physics !Config Key = type_trac !Config Desc = Choix de couplage avec model de chimie INCA ou REPROBUS !Config Def = lmdz !Config Help = !Config 'lmdz' = pas de couplage, pur LMDZ !Config 'lmdz|inca' = model de chime INCA !Config 'lmdz|repr' = model de chime REPROBUS type_trac = 'lmdz' CALL getin('type_trac',type_trac) descrq( 1: 2) = ['LMV','BAK'] descrq(10:20) = ['VL1','VLP','FH1','FH2','VLH',' ','PPM','PPS','PPP',' ','SLP'] descrq(30) = 'PRA' nqo = 0 DO iq=1,nqtot tracers_ico2lmdz(iq)%name = tracs(iq)%name tracers_ico2lmdz(iq)%gen0Name = tracs(iq)%name tracers_ico2lmdz(iq)%phase = tracs(iq)%phase tracers_ico2lmdz(iq)%iadv = tracs(iq)%iadv IF (tracs(iq)%component .eq. "dynamico") then tracers_ico2lmdz(iq)%component='lmdz' ELSE tracers_ico2lmdz(iq)%component=tracs(iq)%component ENDIF if (tracers_ico2lmdz(iq)%iadv .ne. 0 ) tracers_ico2lmdz(iq)%isAdvected=.true. tracers_ico2lmdz(iq)%longName = tracers_ico2lmdz(iq)%name IF(tracers_ico2lmdz(iq)%iadv > 0) tracers_ico2lmdz(iq)%longName=TRIM(tracers_ico2lmdz(iq)%name)//descrq(tracers_ico2lmdz(iq)%iadv) tracers_ico2lmdz(iq)%isInPhysics= delPhase(tracers_ico2lmdz(iq)%gen0Name) /= 'H2O' .OR. tracers_ico2lmdz(iq)%component /= 'lmdz' tracers_ico2lmdz(iq)%iGeneration = 0 ENDDO nqo = COUNT(delPhase(tracs(:)%name) == 'H2O' .AND. tracers_ico2lmdz(:)%component == 'lmdz') !--- Number of water phases isotopes_ico2lmdz(1)%parent='H2O' isotopes_ico2lmdz(1)%phase='gls' isotopes_ico2lmdz(1)%nphas=3 nbtr=nqtot-nqo ALLOCATE(conv_flg(nbtr)) ALLOCATE(pbl_flg(nbtr)) conv_flg(:) = 1 ! convection activated for all tracers pbl_flg(:) = 1 ! boundary layer activated for all tracers CALL init_infotrac_phy ! Initialize physical constant day_length=86400 CALL getin('day_length',day_length) CALL inifis(day_length,radius,g,kappa*cpp,cpp) ! init time annee_ref=2015 CALL getin("anneeref",annee_ref) day_ref=1 CALL getin("dayref",day_ref) physics_timestep=dt*itau_physics run_length=itaumax*dt ndays=NINT(run_length/day_length) day_ini=INT(itau0*dt/day_length)+day_ref start_time= itau0*dt/day_length-INT(itau0*dt/day_length) CALL init_time_lmdz(annee_ref, day_ref, day_ini, start_time, int(ndays), physics_timestep) ! Additional initializations for aquaplanets ! CALL getin("iflag_phys",iflag_phys) ! IF (iflag_phys>=100) THEN ! CALL iniaqua(klon_omp, iflag_phys) ! END IF #ifdef INCA CONTAINS SUBROUTINE init_chem_trac() IMPLICIT NONE CALL Init_chem_inca_trac(nbtr) END SUBROUTINE init_chem_trac SUBROUTINE init_chem_transport() IMPLICIT NONE CALL init_transport(solsym, conv_flg,pbl_flg, hadv_inca, vadv_inca) END SUBROUTINE init_chem_transport #else CONTAINS SUBROUTINE init_chem_trac() IMPLICIT NONE END SUBROUTINE init_chem_trac SUBROUTINE init_chem_transport() IMPLICIT NONE END SUBROUTINE init_chem_transport #endif END SUBROUTINE initialize_physics_omp SUBROUTINE physics USE icosa USE time_mod USE disvert_mod USE transfert_mod USE mpipara USE omp_para USE xios_mod USE wxios USE trace USE distrib_icosa_lmdz_mod, ONLY : transfer_icosa_to_lmdz, transfer_lmdz_to_icosa USE physics_external_mod, ONLY : it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q USE write_field_mod USE checksum_mod USE vorticity_mod USE tracer_icosa_mod USE compute_transport_mod USE, INTRINSIC :: IEEE_ARITHMETIC ! from LMDZ USE mod_phys_lmdz_omp_data, ONLY: klon_omp USE geometry_mod, ONLY : cell_area USE physiq_mod, ONLY: physiq USE icolmdz_param_gravity_wave, ONLY: param_gravity_wave USE mod_phys_lmdz_para, ONLY: reduce_sum, reduce_min USE isotopes_mod, ONLY : iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO, toce, Rdefault, using_iso, initialisation_iso ! ==> LMDZ USE write_field_phy IMPLICIT NONE REAL(rstd),POINTER :: phis(:) REAL(rstd),POINTER :: ps(:) REAL(rstd),POINTER :: theta_rhodz(:,:,:) REAL(rstd),POINTER :: u(:,:) REAL(rstd),POINTER :: wflux(:,:) REAL(rstd),POINTER :: q(:,:,:) REAL(rstd),POINTER :: p(:,:) REAL(rstd),POINTER :: pks(:) REAL(rstd),POINTER :: pk(:,:) REAL(rstd),POINTER :: p_layer(:,:) REAL(rstd),POINTER :: theta(:,:) REAL(rstd),POINTER :: phi(:,:) REAL(rstd),POINTER :: Temp(:,:) REAL(rstd),POINTER :: ulon(:,:) REAL(rstd),POINTER :: ulat(:,:) REAL(rstd),POINTER :: vort(:,:) REAL(rstd),POINTER :: vortc(:,:) REAL(rstd),POINTER :: dulon(:,:) REAL(rstd),POINTER :: dulat(:,:) REAL(rstd),POINTER :: dTemp(:,:) REAL(rstd),POINTER :: dq(:,:,:) REAL(rstd),POINTER :: dps(:) REAL(rstd),POINTER :: duc(:,:,:) INTEGER :: ind,l,k, iq, iq_child REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:) !$OMP THREADPRIVATE(ps_phy) REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:) !$OMP THREADPRIVATE(p_phy) REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:) !$OMP THREADPRIVATE(p_layer_phy) REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:) !$OMP THREADPRIVATE(Temp_phy) REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:) !$OMP THREADPRIVATE(phis_phy) REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:) !$OMP THREADPRIVATE(phi_phy) REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:) !$OMP THREADPRIVATE(ulon_phy) REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:) !$OMP THREADPRIVATE(ulat_phy) REAL(rstd),ALLOCATABLE,SAVE :: rot_phy(:,:) !$OMP THREADPRIVATE(rot_phy) REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:) !$OMP THREADPRIVATE(q_phy) REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:) !$OMP THREADPRIVATE(wflux_phy) REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:) !$OMP THREADPRIVATE(dulon_phy) REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:) !$OMP THREADPRIVATE(dulat_phy) REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:) !$OMP THREADPRIVATE(dTemp_phy) REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:) !$OMP THREADPRIVATE(dq_phy) REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:) !$OMP THREADPRIVATE(dps_phy) REAL(rstd) :: dtphy LOGICAL :: debut LOGICAL :: lafin LOGICAL,SAVE :: first=.TRUE. !$OMP THREADPRIVATE(first) INTEGER, SAVE :: iter !$OMP THREADPRIVATE(iter) REAL(rstd),ALLOCATABLE,SAVE :: q_phy_old(:,:,:) REAL :: qmin_omp,qmin_glo REAL :: q_ave_bis(llm,nqtot) IF(first) THEN debut=.TRUE. iter = 1 ELSE debut=.FALSE. iter=iter+1 ENDIF IF(it-itau0>=itaumax) THEN lafin=.TRUE. ELSE lafin=.FALSE. ENDIF IF (first) THEN first=.FALSE. CALL init_message(f_u,req_e1_vect,req_u) CALL init_message(f_vort,req_z1_scal,req_z) ALLOCATE(ps_phy(klon_omp)) ALLOCATE(p_phy(klon_omp,llm+1)) ALLOCATE(p_layer_phy(klon_omp,llm)) ALLOCATE(Temp_phy(klon_omp,llm)) ALLOCATE(phis_phy(klon_omp)) ALLOCATE(phi_phy(klon_omp,llm)) ALLOCATE(ulon_phy(klon_omp,llm)) ALLOCATE(ulat_phy(klon_omp,llm)) ALLOCATE(rot_phy(klon_omp,llm)) ALLOCATE(q_phy(klon_omp,llm,nqtot)) ALLOCATE(wflux_phy(klon_omp,llm)) ALLOCATE(dulon_phy(klon_omp,llm)) ALLOCATE(dulat_phy(klon_omp,llm)) ALLOCATE(dTemp_phy(klon_omp,llm)) ALLOCATE(dq_phy(klon_omp,llm,nqtot)) ALLOCATE(dps_phy(klon_omp)) ALLOCATE(q_phy_old(klon_omp,llm,nqtot)) !$OMP BARRIER ! CALL param_gravity_wave ENDIF !$OMP MASTER ! CALL update_calendar(it) !$OMP END MASTER !$OMP BARRIER dtphy=itau_physics*dt CALL transfert_message(f_u,req_u) DO ind=1,ndomain IF (assigned_domain(ind)) THEN CALL swap_dimensions(ind) CALL swap_geometry(ind) u=f_u(ind) vort=f_vort(ind) CALL compute_vorticity(u,vort) ENDIF ENDDO CALL transfert_message(f_vort,req_z) DO ind=1,ndomain CALL swap_dimensions(ind) IF (assigned_domain(ind)) THEN CALL swap_geometry(ind) phis=f_phis(ind) ps=f_ps(ind) theta_rhodz=f_theta_rhodz(ind) u=f_u(ind) q=f_q(ind) wflux=f_wflux(ind) p=f_p(ind) pks=f_pks(ind) pk=f_pk(ind) p_layer=f_p_layer(ind) theta=f_theta(ind) phi=f_phi(ind) Temp=f_Temp(ind) ulon=f_ulon(ind) ulat=f_ulat(ind) vort=f_vort(ind) vortc=f_vortc(ind) CALL grid_icosa_to_physics ENDIF ENDDO !$OMP BARRIER !$OMP MASTER CALL SYSTEM_CLOCK(start_clock) !$OMP END MASTER CALL trace_start("physic") ! CALL trace_off() ! CALL writeField("p_in",f_p) ! CALL writeField("p_layer_in",f_p_layer) ! CALL writeField("phi_in",f_phi) ! CALL writeField("phis_in",f_phis) ! CALL writeField("ulon_in",f_ulon) ! CALL writeField("ulat_in",f_ulat) ! CALL writeField("Temp_in",f_Temp) ! CALL writeField("q_in",f_q) ! CALL writeField("wflux_in",f_wflux) ! CALL writeField("vortc",f_vortc) ! CALL checksum(f_p) ! CALL checksum(f_p_layer) ! CALL checksum(f_phi) ! CALL checksum(f_phis) ! CALL checksum(f_ulon) ! CALL checksum(f_ulat) ! CALL checksum(f_Temp) ! CALL checksum(f_q) ! CALL checksum(f_wflux) CALL transfer_icosa_to_lmdz(f_p , p_phy) CALL transfer_icosa_to_lmdz(f_p_layer, p_layer_phy) CALL transfer_icosa_to_lmdz(f_phi , phi_phy) CALL transfer_icosa_to_lmdz(f_phis , phis_phy ) CALL transfer_icosa_to_lmdz(f_ulon , ulon_phy ) CALL transfer_icosa_to_lmdz(f_ulat , ulat_phy) CALL transfer_icosa_to_lmdz(f_vortc , rot_phy) CALL transfer_icosa_to_lmdz(f_Temp , Temp_phy) CALL transfer_icosa_to_lmdz(f_q , q_phy) CALL transfer_icosa_to_lmdz(f_wflux , wflux_phy) DO l=1,llm wflux_phy(:,l) = - wflux_phy(:,l)*cell_area(:) phi_phy(:,l)=phi_phy(:,l)-phis_phy(:) ENDDO DO iq=1, nqtot IF (.NOT. tracers(iq)%has_parent) THEN DO iq_child=1,tracers(iq)%nb_children CALL tracers_dyn2phys(q_phy, tracers(iq)%children(iq_child), iq) ENDDO ENDIF ENDDO CALL wxios_set_context() ! Ehouarn: rot_phy() not implemented!! Set it to zero for now ! rot_phy(:,:)=0 CALL physiq(klon_omp, llm, debut, lafin, dtphy, & p_phy, p_layer_phy, phi_phy, phis_phy, presnivs, & ulon_phy, ulat_phy, rot_phy, Temp_phy, q_phy, wflux_phy, & dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy) CALL compute_q_ave(q_phy, q_ave) DO iq=1, nqtot q_phy(:,:,iq) = q_phy(:,:,iq) + dtphy*dq_phy(:,:,iq) ENDDO ! IF (is_master) PRINT*, "==================================================" ! IF (is_master) PRINT*, "==================================================" ! CALL compute_q_ave(q_phy, q_ave_bis) ! DO iq=1, nqtot ! qmin_omp = MINVAL(q_phy(:,:,iq)) ! CALL reduce_min(qmin_omp, qmin_glo) ! IF (is_master) PRINT*, "after 0 : MINVAL of ",TRIM(tracers(iq)%name), qmin_glo ! IF (is_master) PRINT*, "AVERAGE of ",TRIM(tracers(iq)%name), q_ave_bis(:,iq) ! ENDDO CALL compute_qminimum(q_phy, q_ave, p_phy) ! DO iq=1, nqtot ! IF (.NOT. tracers(iq)%has_parent) THEN ! DO iq_child=1,tracers(iq)%nb_children ! CALL tracers_normalize(q_phy, tracers(iq)%children(iq_child), iq) ! ENDDO ! DO l=1,llm ! DO k=1,klon_omp ! IF ( q_phy(k,l,iq) < q_epsilon*q_ave(l,iq)) q_phy(k,l,iq) = q_epsilon*q_ave(l,iq) ! ENDDO ! ENDDO ! ENDIF ! ENDDO ! IF (is_master) PRINT*, "==================================================" ! IF (is_master) PRINT*, "==================================================" ! ! CALL compute_q_ave(q_phy, q_ave) ! ! DO iq=1, nqtot ! qmin_omp = MINVAL(q_phy(:,:,iq)) ! CALL reduce_min(qmin_omp, qmin_glo) ! IF (is_master) PRINT*, "after 1 MINVAL of ",TRIM(tracers(iq)%name), qmin_glo ! IF (is_master) PRINT*, "AVERAGE of ",TRIM(tracers(iq)%name), q_ave(:,iq) ! IF (ANY(IEEE_IS_NAN(q_ave(:,iq)))) THEN ! CALL WriteField_phy("debug_"//TRIM(tracers(iq)%name),q_phy(:,:,iq),llm) ! IF (is_master) PRINT *,'NAN DETECTED' ! STOP ! ENDIF ! ENDDO DO iq=1, nqtot IF (.NOT. tracers(iq)%has_parent) THEN DO iq_child=1,tracers(iq)%nb_children CALL tracers_phys2dyn(q_phy, tracers(iq)%children(iq_child), iq) ENDDO ENDIF ENDDO CALL transfer_lmdz_to_icosa(dulon_phy, f_dulon ) CALL transfer_lmdz_to_icosa(dulat_phy, f_dulat ) CALL transfer_lmdz_to_icosa(dTemp_phy, f_dTemp ) CALL transfer_lmdz_to_icosa(q_phy , f_q ) CALL transfer_lmdz_to_icosa(dps_phy , f_dps ) ! CALL writeField("dulon_out",f_dulon) ! CALL writeField("dulat_out",f_dulat) ! CALL writeField("dTemp_out",f_dTemp) ! CALL writeField("dq_out",f_dq) ! CALL writeField("dps_out",f_dps) ! CALL checksum(f_dulon) ! CALL checksum(f_dulat) ! CALL checksum(f_dTemp) ! CALL checksum(f_dq) ! CALL checksum(f_dps) CALL send_message(f_dps,req_dps0) CALL send_message(f_dulon,req_dulon0) CALL send_message(f_dulat,req_dulat0) CALL send_message(f_dTemp,req_dTemp0) CALL send_message(f_q,req_dq0) CALL wait_message(req_dps0) CALL wait_message(req_dulon0) CALL wait_message(req_dulat0) CALL wait_message(req_dTemp0) CALL wait_message(req_dq0) ! CALL trace_on() CALL trace_end("physic") !$OMP MASTER CALL SYSTEM_CLOCK(stop_clock) count_clock=count_clock+stop_clock-start_clock !$OMP END MASTER !$OMP BARRIER DO ind=1,ndomain CALL swap_dimensions(ind) IF (assigned_domain(ind)) THEN CALL swap_geometry(ind) theta_rhodz=f_theta_rhodz(ind) u=f_u(ind) q=f_q(ind) ps=f_ps(ind) dulon=f_dulon(ind) dulat=f_dulat(ind) Temp=f_temp(ind) dTemp=f_dTemp(ind) dq=f_dq(ind) dps=f_dps(ind) duc=f_duc(ind) p=f_p(ind) pks=f_pks(ind) pk=f_pk(ind) CALL grid_physics_to_icosa ENDIF ENDDO !$OMP BARRIER CALL xios_set_context CONTAINS SUBROUTINE grid_icosa_to_physics USE pression_mod USE exner_mod USE theta2theta_rhodz_mod USE geopotential_mod USE wind_from_lonlat_mod USE omp_para IMPLICIT NONE REAL(rstd) :: uc(3) INTEGER :: i,j,ij,l ! compute pression DO l = ll_begin,ll_endp1 DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i p(ij,l) = ap(l) + bp(l) * ps(ij) ENDDO ENDDO ENDDO !$OMP BARRIER ! compute exner IF (is_omp_first_level) THEN DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i pks(ij) = cpp * ( ps(ij)/preff ) ** kappa ENDDO ENDDO ENDIF ! 3D : pk DO l = ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa ENDDO ENDDO ENDDO !$OMP BARRIER ! compute theta, temperature and pression at layer DO l = ll_begin, ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g) Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp p_layer(ij,l)=preff*(pk(ij,l)/cpp)**(1./kappa) ENDDO ENDDO ENDDO !!! Compute geopotential ! for first layer IF (is_omp_first_level) THEN DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) ENDDO ENDDO ENDIF !!-> implicit flush on phi(:,1) !$OMP BARRIER ! for other layers DO l = ll_beginp1, ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i phi(ij,l) = 0.5 * ( theta(ij,l) + theta(ij,l-1) ) & * ( pk(ij,l-1) - pk(ij,l) ) ENDDO ENDDO ENDDO !$OMP BARRIER IF (is_omp_first_level) THEN DO l = 2, llm DO j=jj_begin,jj_end ! ---> Bug compilo intel ici en openmp ! ---> Couper la boucle IF (j==jj_end+1) PRINT*,"this message must not be printed" DO i=ii_begin,ii_end ij=(j-1)*iim+i phi(ij,l) = phi(ij,l)+ phi(ij,l-1) ENDDO ENDDO ENDDO ! --> IMPLICIT FLUSH on phi --> non ENDIF ! compute wind centered lon lat compound DO l=ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i uc(:)=1/Ai(ij)* & ( ne(ij,right)*u(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:)) & + ne(ij,rup)*u(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:)) & + ne(ij,lup)*u(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:)) & + ne(ij,left)*u(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:)) & + ne(ij,ldown)*u(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))& + ne(ij,rdown)*u(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:))) ulon(ij,l)=sum(uc(:)*elon_i(ij,:)) ulat(ij,l)=sum(uc(:)*elat_i(ij,:)) ENDDO ENDDO ENDDO ! compute centered vorticity DO l=ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i vortc(ij,l) = Riv(ij,vup) * vort(ij+z_up,l) + & Riv(ij,vlup) * vort(ij+z_lup,l) + & Riv(ij,vldown)* vort(ij+z_ldown,l) + & Riv(ij,vdown) * vort(ij+z_down,l) + & Riv(ij,vrdown)* vort(ij+z_rdown,l) + & Riv(ij,vrup) * vort(ij+z_rup,l) ENDDO ENDDO ENDDO !$OMP BARRIER END SUBROUTINE grid_icosa_to_physics SUBROUTINE grid_physics_to_icosa USE theta2theta_rhodz_mod USE omp_para IMPLICIT NONE INTEGER :: i,j,ij,l,iq DO l=ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:) ENDDO ENDDO ENDDO DO l=ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) ) u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) ) u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) ) ENDDO ENDDO ENDDO DO l=ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l) ENDDO ENDDO ENDDO ! DO iq=1,nqtot ! DO l=ll_begin,ll_end ! DO j=jj_begin,jj_end ! DO i=ii_begin,ii_end ! ij=(j-1)*iim+i ! q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq) ! ENDDO ! ENDDO ! ENDDO ! ENDDO !$OMP BARRIER IF (is_omp_first_level) THEN DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i ps(ij)=ps(ij)+ dtphy * dps(ij) ENDDO ENDDO ENDIF ! CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0) ! compute pression !$OMP BARRIER DO l = ll_begin,ll_endp1 DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i p(ij,l) = ap(l) + bp(l) * ps(ij) ENDDO ENDDO ENDDO !$OMP BARRIER ! compute exner IF (is_omp_first_level) THEN DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i pks(ij) = cpp * ( ps(ij)/preff ) ** kappa ENDDO ENDDO ENDIF ! 3D : pk DO l = ll_begin,ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa ENDDO ENDDO ENDDO !$OMP BARRIER ! compute theta, temperature and pression at layer DO l = ll_begin, ll_end DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i theta_rhodz(ij,l,1) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp ) ENDDO ENDDO ENDDO END SUBROUTINE grid_physics_to_icosa RECURSIVE SUBROUTINE tracers_dyn2phys(q, iq, iq_parent) USE tracer_icosa_mod IMPLICIT NONE REAL(rstd),INTENT(INOUT) :: q(:,:,:) INTEGER,INTENT(IN) :: iq, iq_parent INTEGER :: iq_child q(:,:,iq) = q(:,:,iq)*q(:,:,iq_parent) DO iq_child=1,tracers(iq)%nb_children CALL tracers_dyn2phys(q, tracers(iq)%children(iq_child), iq) ENDDO END SUBROUTINE tracers_dyn2phys RECURSIVE SUBROUTINE tracers_phys2dyn(q, iq, iq_parent) USE tracer_icosa_mod USE compute_transport_mod USE mod_phys_lmdz_para USE mod_phys_lmdz_omp_data, ONLY: klon_omp USE dimphy USE mod_grid_phy_lmdz IMPLICIT NONE REAL(rstd),INTENT(INOUT) :: q(klon_omp,nbp_lev,nqtot) INTEGER,INTENT(IN) :: iq, iq_parent INTEGER :: iq_child INTEGER :: lij DO iq_child=1,tracers(iq)%nb_children CALL tracers_phys2dyn(q, tracers(iq)%children(iq_child), iq) ENDDO q(:,:,iq) = q(:,:,iq) / q(:,:,iq_parent) END SUBROUTINE tracers_phys2dyn RECURSIVE SUBROUTINE tracers_normalize(q, iq, iq_parent) USE tracer_icosa_mod USE compute_transport_mod USE mod_phys_lmdz_para USE mod_phys_lmdz_omp_data, ONLY: klon_omp USE dimphy USE mod_grid_phy_lmdz IMPLICIT NONE REAL(rstd),INTENT(INOUT) :: q(klon_omp,nbp_lev,nqtot) INTEGER,INTENT(IN) :: iq, iq_parent INTEGER :: iq_child INTEGER :: l, k DO iq_child=1,tracers(iq)%nb_children CALL tracers_normalize(q, tracers(iq)%children(iq_child), iq) ENDDO DO l=1,klev DO k=1,klon IF ( q(k,l,iq_parent) >= q_epsilon*q_ave(l, iq_parent) ) THEN IF ( q(k,l,iq) < q_epsilon*q_ave(l, iq_parent) ) THEN q(k,l,iq) = q_ave(l, iq)/q_ave(l, iq_parent) * q(k,l,iq_parent) ENDIF ELSE q(k,l,iq) = q_ave(l, iq) * q_epsilon ENDIF ENDDO ENDDO END SUBROUTINE tracers_normalize SUBROUTINE compute_q_ave(q, q_ave) USE dimphy USE tracer_icosa_mod USE mod_phys_lmdz_para USE mod_grid_phy_lmdz IMPLICIT NONE REAL(rstd),INTENT(IN) :: q(klon, klev ,nqtot) REAL(rstd),INTENT(OUT) :: q_ave(klev, nqtot) REAL(rstd) :: q_ave_omp(klev, nqtot) INTEGER :: iq, l , k DO iq=1,nqtot DO l=1,klev q_ave_omp(l,iq) = 0. DO k=1,klon q_ave_omp(l,iq) = q_ave_omp(l,iq) + q(k,l,iq) ENDDO ENDDO ENDDO CALL reduce_sum(q_ave_omp, q_ave) CALL bcast(q_ave) q_ave(:,:) = q_ave(:,:) / klon_glo WHERE (q_ave < q_minimum) q_ave = q_minimum END SUBROUTINE compute_q_ave SUBROUTINE compute_qminimum(q, q_ave, p ) USE dimphy USE tracer_icosa_mod IMPLICIT NONE REAL(rstd),INTENT(INOUT) :: q(klon, klev ,nqtot) REAL(rstd),INTENT(IN) :: q_ave(klev, nqtot) REAL(rstd),INTENT(IN) :: p(klon, klev) REAL(rstd) :: deltap(klon, klev) INTEGER :: iq, iq_g INTEGER :: iq_H2O_g, iq_H2O_l, iq_H216O_g, iq_H216O_l , iq_H218O_g, iq_H218O_l, iq_HDO_g, iq_HDO_l DO iq=1, nqtot SELECT CASE (TRIM(tracers(iq)%name)) CASE ("H2O_g") ; iq_H2O_g = iq CASE ("H2O_l") ; iq_H2O_l = iq CASE ("H216O_g") ; iq_H216O_g = iq CASE ("H216O_l") ; iq_H216O_l = iq CASE ("H218O_g") ; iq_H218O_g = iq CASE ("H218O_l") ; iq_H218O_l = iq CASE ("HDO_g") ; iq_HDO_g = iq CASE ("HDO_l") ; iq_HDO_l = iq CASE DEFAULT ; CYCLE END SELECT ENDDO DO iq=1, nqtot IF (iq==iq_H2O_l) THEN iq_g = iq_H2O_g ELSE IF (iq==iq_H216O_l) THEN iq_g = iq_H216O_g ELSE IF (iq==iq_H218O_l) THEN iq_g = iq_H218O_g ELSE IF (iq==iq_HDO_l) THEN iq_g = iq_HDO_g ELSE CYCLE ENDIF DO l=1,klev DO k=1,klon IF ( q(k,l,iq) <= q_epsilon*q_ave(l,iq)) THEN q(k,l,iq_g) = q(k,l,iq_g) - (q_epsilon*q_ave(l,iq) - q(k,l,iq)) q(k,l,iq) = q_epsilon*q_ave(l,iq) ENDIF ENDDO ENDDO ENDDO DO l=1,klev DO k=1,klon deltap(k,l) = p(k,l)-p(k,l+1) ENDDO ENDDO DO iq=1, nqtot DO l=klev,2 DO k=1,klon IF ( q(k,l,iq) < q_epsilon*q_ave(l,iq) ) THEN q(k,l-1,iq) = q(k,l-1,iq) - (q_epsilon*q_ave(l,iq) - q(k,l,iq))*deltap(k,l)/deltap(k,l-1) q(k,l,iq) = q_epsilon*q_ave(l,iq) ENDIF ENDDO ENDDO ENDDO DO iq=1, nqtot DO l=1,klev-1 DO k=1,klon IF ( q(k,l,iq) < q_epsilon*q_ave(l,iq) ) THEN q(k,l+1,iq) = q(k,l+1,iq) - (q_epsilon*q_ave(l,iq) - q(k,l,iq))*deltap(k,l)/deltap(k,l+1) q(k,l,iq) = q_epsilon*q_ave(l,iq) ENDIF ENDDO DO k=1,klon IF (q(k,klev,iq) < q_epsilon*q_ave(klev,iq)) q(k,klev,iq) = q_epsilon*q_ave(klev,iq) ENDDO ENDDO ENDDO END SUBROUTINE compute_qminimum END SUBROUTINE physics END MODULE interface_icosa_lmdz_mod