Index: trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3948) +++ trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3949) @@ -204,5 +204,6 @@ !! Microphysics-specific variables logical,save :: callmufi, call_haze_prod_pCH4 -!$OMP THREADPRIVATE(callmufi,call_haze_prod_pCH4) + logical,save :: callmuclouds +!$OMP THREADPRIVATE(callmufi,call_haze_prod_pCH4,callmuclouds) real,save :: haze_p_prod, haze_tx_prod, haze_rc_prod real,save :: haze_rm, haze_df, haze_rho Index: trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90 (revision 3948) +++ trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90 (revision 3949) @@ -18,5 +18,6 @@ ! Set in inifis_mod character(LEN=100),save :: config_mufi ='datagcm/microphysics/config.cfg' -!$OMP THREADPRIVATE(config_mufi) + character(LEN=100),save :: mugasflux_file ='datagcm/microphysics/Species_Pluto_PLavvas_Flux.txt' +!$OMP THREADPRIVATE(config_mufi,mugasflux_file) character(len=300),save :: aersprop_file character(len=300),save :: aerfprop_file Index: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3948) +++ trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3949) @@ -14,5 +14,5 @@ use radii_mod, only: radfixed use datafile_mod, only: datadir,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file, & - config_mufi, aersprop_file, aerfprop_file + config_mufi, mugasflux_file, aersprop_file, aerfprop_file use comdiurn_h, only: sinlat, coslat, sinlon, coslon use comgeomfi_h, only: totarea, totarea_planet @@ -764,4 +764,9 @@ if (is_master) write(*,*)" callmufi = ",callmufi + if (is_master) write(*,*) "Run with or without microphysical clouds?" + callmuclouds=.false. ! default value + call getin_p("callmuclouds",callmuclouds) + if (is_master) write(*,*)" callmuclouds = ",callmuclouds + ! sanity check if (callmufi.and.(.not.tracer)) then @@ -770,9 +775,19 @@ stop endif + if (callmuclouds.and.(.not.callmufi)) then + print*,"You are running microphysical clouds without microphysics" + print*,"Please start again with callmufi =.true." + stop + endif if (is_master) write(*,*) "Path to microphysical config file?" config_mufi='datagcm/microphysics/config.cfg' ! default value call getin_p("config_mufi",config_mufi) - if (is_master) write(*,*)" config_mufi = ",config_mufi + if (is_master) write(*,*) trim(rname)//" config_mufi = ",config_mufi + + if (is_master) write(*,*) "Condensable gas fluxes datafile" + mugasflux_file='None' ! default file + call getin_p("mugasflux_file",mugasflux_file) + if (is_master) write(*,*) trim(rname)//" mugasflux_file = ",trim(mugasflux_file) if (is_master) write(*,*) "Spherical aerosol optical properties datafile" @@ -1439,4 +1454,7 @@ call abort_physic(rname, 'if microphysics is on, naerkind must be > 1!', 1) endif + if ((callmufi).and.(callmuclouds).and..not.(naerkind.gt.2)) then + call abort_physic(rname, 'if microphysical clouds are on, naerkind must be > 2!', 1) + endif if (.not.(callmufi.or.haze).and.(optichaze)) then call abort_physic(rname, 'if microphysics and haze are off, optichaze must be deactivated!', 1) Index: trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90 (revision 3948) +++ trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90 (revision 3949) @@ -3,5 +3,6 @@ use surfdat_h, ONLY: dryness USE tracer_h - USE callkeys_mod, only: optichaze,nb_monomer,haze,fractal,fasthaze,rad_haze,callmufi + USE callkeys_mod, only: optichaze,nb_monomer,haze,fractal,fasthaze,rad_haze,& + callmufi,callmuclouds USE recombin_corrk_mod, ONLY: ini_recombin USE mod_phys_lmdz_para, only: is_master, bcast @@ -150,7 +151,12 @@ igcm_ch4_gas=0 igcm_ch4_ice=0 - igcm_prec_haze=0 igcm_co_gas=0 igcm_co_ice=0 + igcm_C2H2_mugas=0 + igcm_C2H6_mugas=0 + igcm_C4H2_mugas=0 + igcm_C6H6_mugas=0 + igcm_HCN_mugas=0 + igcm_prec_haze=0 nqhaze(:)=0 @@ -225,4 +231,36 @@ write(*,*) 'Tracer ',count,' = co ice' endif +! Microphysics related tracers + if (noms(iq).eq."C2H2_mugas") then + igcm_C2H2_mugas=iq + mmol(igcm_C2H2_mugas)=26.04 + count=count+1 + write(*,*) 'Tracer ',count,' = C2H2 mugas' + endif + if (noms(iq).eq."C2H6_mugas") then + igcm_C2H6_mugas=iq + mmol(igcm_C2H6_mugas)=30.07 + count=count+1 + write(*,*) 'Tracer ',count,' = C2H6 mugas' + endif + if (noms(iq).eq."C4H2_mugas") then + igcm_C4H2_mugas=iq + mmol(igcm_C4H2_mugas)=50.05 + count=count+1 + write(*,*) 'Tracer ',count,' = C4H2 mugas' + endif + if (noms(iq).eq."C6H6_mugas") then + igcm_C6H6_mugas=iq + mmol(igcm_C6H6_mugas)=78.11 + count=count+1 + write(*,*) 'Tracer ',count,' = C6H6 mugas' + endif + if (noms(iq).eq."HCN_mugas") then + igcm_HCN_mugas=iq + mmol(igcm_HCN_mugas)=27.03 + count=count+1 + write(*,*) 'Tracer ',count,' = HCN mugas' + endif +! Haze tracers if (noms(iq).eq."prec_haze") then igcm_prec_haze=iq @@ -293,31 +331,5 @@ enddo ! of do iq=1,nq - ! ! 3. find condensable traceurs different from h2o and n2 - ! do iq=1,nq - ! if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then - ! count=count+1 - ! endif - ! if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then - ! count=count+1 - ! endif - - ! enddo ! of do iq=1,nq - - ! check that we identified all tracers: - if (count.ne.nq) then - write(*,*) "initracer: found only ",count," tracers" - write(*,*) " expected ",nq - do iq=1,count - write(*,*)' ',iq,' ',trim(noms(iq)) - enddo -! stop - else - write(*,*) "initracer: found all expected tracers, namely:" - do iq=1,nq - write(*,*)' ',iq,' ',trim(noms(iq)) - enddo - endif - - ! Compute number of microphysics tracers: + ! 3. Find microphysics tracers: ! By convention they all have the prefix "mu_" (case sensitive !) nmicro = 0 @@ -331,16 +343,26 @@ ENDDO - ! Checking the expected number of tracers: - ! Microphysics moment model: nmicro = 4 - IF (nmicro < 4) THEN + ! Checking the expected number of microphysical tracers: + ! No cloud: nmicro = 4 aer. + ! Clouds: nmicro = 4 aer. + 2 ccn + 1(+) ices + IF (callmuclouds) THEN + IF (nmicro < 7) THEN WRITE(*,*) "initracer:error:"," Inconsistent number of microphysical tracers" - WRITE(*,*) "expected at least 4 tracers,", nmicro, " given" + WRITE(*,*) "expected at least 7 tracers (clouds: on),", nmicro, " given" CALL abort - ELSE IF (nmicro > 4) THEN + ENDIF + + ELSE + IF (nmicro < 4) THEN + WRITE(*,*) "initracer:error:"," Inconsistent number of microphysical tracers" + WRITE(*,*) "expected at least 4 tracers (clouds: off),", nmicro, " given" + CALL abort + ELSE IF (nmicro > 4) THEN WRITE(*,*) "!!! WARNING !!! initracer: I was expecting only four tracers, you gave me more." CALL abort - ENDIF - - ! microphysics indexes share the same values than original tracname. + ENDIF + ENDIF ! end of callmuclouds + + ! Microphysics indexes share the same values than original tracname. IF (.NOT.ALLOCATED(micro_indx)) ALLOCATE(micro_indx(nmicro)) j = 1 @@ -358,4 +380,18 @@ ENDIF ! end of callmufi + ! Check that we identified all tracers: + if (count.ne.nq) then + write(*,*) "initracer: found only ",count," tracers" + write(*,*) " expected ",nq + do iq=1,count + write(*,*)' ',iq,' ',trim(noms(iq)) + enddo + else + write(*,*) "initracer: found all expected tracers, namely:" + do iq=1,nq + write(*,*)' ',iq,' ',trim(noms(iq)) + enddo + endif + ! Get data of tracers. Need to rewind traceur.def first if (is_master) then @@ -376,10 +412,4 @@ if (is_master) close(407) - ! Calculate number of microphysical tracer - write(*,*) 'Number of microphysical tracer nmicro = ',nmicro - IF (callmufi) THEN - call dumptracers(micro_indx) - ENDIF - ! Processing modern traceur options if(moderntracdef) then @@ -401,7 +431,9 @@ if (callmufi) then if (optichaze) then - iaero_haze = 2 - write(*,*) 'Microphysical moment model' - write(*,*) '--- number of haze aerosol = ', iaero_haze + if (callmuclouds) then + iaero_haze = 3 + else + iaero_haze = 2 + endif ! end of callmuclouds endif ! end optichaze @@ -431,5 +463,5 @@ (first one in traceur.def)" endif - enddo + enddo endif ! end optichaze endif ! end callmufi or haze Index: trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90 (revision 3948) +++ trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90 (revision 3949) @@ -64,4 +64,9 @@ haze_rm,haze_df,haze_rho,rad,g,air_rad,mugaz*1e-3, & config_mufi) - + + ! Sanity check for haze model initialization: + ! ------------------------------------------- + write(*,*) 'Number of microphysical tracer nmicro = ',nmicro + call dumptracers(micro_indx) + end subroutine inimufi Index: trunk/LMDZ.PLUTO/libf/phypluto/mugas_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/mugas_mod.F90 (revision 3949) +++ trunk/LMDZ.PLUTO/libf/phypluto/mugas_mod.F90 (revision 3949) @@ -0,0 +1,175 @@ +subroutine mugas_prof(ngrid,nlayer,nq,zzlay,zzlev,pplay,pt,pdqchem) + use datafile_mod + use comcstfi_mod, only: r + use tracer_h, only: noms,igcm_C2H2_mugas,igcm_C2H6_mugas,& + igcm_C4H2_mugas,igcm_C6H6_mugas,igcm_HCN_mugas + use mod_phys_lmdz_para, only : is_master + + implicit none + +!================================================================== +! Purpose +! ------- +! Get fixed fluxes of microphysical condensable gas +! from txt file. +! +! Inputs +! ------ +! ngrid | Number of atmospheric columns. +! nlayer | Number of atmospheric layers. +! nq | Number of tracers. +! zzlay | Altitude at the middle of the atmospheric layers. +! zzlev | Altitude at the atmospheric layer boundaries. +! pplay | Mid-layer pressure (Pa) +! pt | Temperature (K). +! +! Outputs +! ------- +! +! Both +! ---- +! pdqchem ! Tracer tendencies for condensable gases through muphi (kg/kg_of_air/s). +! +! Authors +! ------- +! Bruno de Batz de Trenquelléon (2025) +!================================================================== + + ! Arguments: + !~~~~~~~~~~~ + integer,intent(in) :: ngrid + integer,intent(in) :: nlayer + integer,intent(in) :: nq + real,intent(in) :: zzlay(ngrid,nlayer) + real,intent(in) :: zzlev(ngrid,nlayer+1) + real,intent(in) :: pplay(ngrid,nlayer) + real,intent(in) :: pt(ngrid,nlayer) + + real,intent(inout) :: pdqchem(ngrid,nlayer,nq) + + ! Local variables: + !~~~~~~~~~~~~~~~~~ + integer :: ig, iq, l, ifine + + real :: dz(ngrid,nlayer) + real :: zdq_HCN(ngrid,nlayer) + real :: zdq_C4H2(ngrid,nlayer) + real :: zdq_C6H6(ngrid,nlayer) + real :: zdq_C2H2(ngrid,nlayer) + real :: zdq_C2H6(ngrid,nlayer) + + integer Nfine + parameter(Nfine=999) + character(len=100) :: file_path + character(len=100) :: file_name + real,save,allocatable :: Altdata(:) + real,save,allocatable :: Pdata(:) + real,save,allocatable :: Tdata(:) + real,save,allocatable :: C2H2_fluxdata(:) + real,save,allocatable :: C2H6_fluxdata(:) + real,save,allocatable :: C4H2_fluxdata(:) + real,save,allocatable :: C6H6_fluxdata(:) + real,save,allocatable :: HCN_fluxdata(:) + + LOGICAL,SAVE :: firstcall=.true. +!$OMP THREADPRIVATE(firstcall) + + ! 0. Initialization: + !~~~~~~~~~~~~~~~~~~~ + dz(:,:) = zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer) + + zdq_C2H2(:,:) = 0.d0 + zdq_C2H6(:,:) = 0.d0 + zdq_C4H2(:,:) = 0.d0 + zdq_C6H6(:,:) = 0.d0 + zdq_HCN(:,:) = 0.d0 + + ! 1. Read data: + !~~~~~~~~~~~~~~ + if (firstcall) then + firstcall=.false. + + if (mugasflux_file/='None') then + file_name = mugasflux_file + else + STOP "No filename given for condensable gas fluxes profile. Need to set mugasflux_file..." + endif +!$OMP MASTER + + if(.not.allocated(Altdata)) allocate(Altdata(Nfine)) + if(.not.allocated(Pdata)) allocate(Pdata(Nfine)) + if(.not.allocated(Tdata)) allocate(Tdata(Nfine)) + if(.not.allocated(C2H2_fluxdata)) allocate(C2H2_fluxdata(Nfine)) + if(.not.allocated(C2H6_fluxdata)) allocate(C2H6_fluxdata(Nfine)) + if(.not.allocated(C4H2_fluxdata)) allocate(C4H2_fluxdata(Nfine)) + if(.not.allocated(C6H6_fluxdata)) allocate(C6H6_fluxdata(Nfine)) + if(.not.allocated(HCN_fluxdata)) allocate(HCN_fluxdata(Nfine)) + + file_path=file_name + open(224,file=file_path,form='formatted') + do ifine = 1, Nfine + read(224,*) Altdata(ifine), Pdata(ifine), Tdata(ifine), HCN_fluxdata(ifine), & + C4H2_fluxdata(ifine), C6H6_fluxdata(ifine), C2H2_fluxdata(ifine), C2H6_fluxdata(ifine) + enddo + close(224) +!$OMP END MASTER +!$OMP BARRIER + endif + + ! 2. Interpolate on the model vertical grid: + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + do iq = 1, nq + + ! C2H2 + if (noms(iq).eq."C2H2_mugas") then + do ig = 1, ngrid + call interp_line(Altdata*1000.,C2H2_fluxdata,Nfine,zzlay(ig,:),zdq_C2H2(ig,:),nlayer) + do l = 1, nlayer + pdqchem(ig,l,igcm_C2H2_mugas) = zdq_C2H2(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) + enddo + enddo + endif + + ! C2H6 + if (noms(iq).eq."C2H6_mugas") then + do ig = 1, ngrid + call interp_line(Altdata*1000.,C2H6_fluxdata,Nfine,zzlay(ig,:),zdq_C2H6(ig,:),nlayer) + do l = 1, nlayer + pdqchem(ig,l,igcm_C2H6_mugas) = zdq_C2H6(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) + enddo + enddo + endif + + ! C4H2 + if (noms(iq).eq."C4H2_mugas") then + do ig = 1, ngrid + call interp_line(Altdata*1000.,C4H2_fluxdata,Nfine,zzlay(ig,:),zdq_C4H2(ig,:),nlayer) + do l = 1, nlayer + pdqchem(ig,l,igcm_C4H2_mugas) = zdq_C4H2(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) + enddo + enddo + endif + + ! C6H6 + if (noms(iq).eq."C6H6_mugas") then + do ig = 1, ngrid + call interp_line(Altdata*1000.,C6H6_fluxdata,Nfine,zzlay(ig,:),zdq_C6H6(ig,:),nlayer) + do l = 1, nlayer + pdqchem(ig,l,igcm_C6H6_mugas) = zdq_C6H6(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) + enddo + enddo + endif + + ! HCN + if (noms(iq).eq."HCN_mugas") then + do ig = 1, ngrid + call interp_line(Altdata*1000.,HCN_fluxdata,Nfine,zzlay(ig,:),zdq_HCN(ig,:),nlayer) + do l = 1, nlayer + pdqchem(ig,l,igcm_HCN_mugas) = zdq_HCN(ig,l) / dz(ig,l) / (pplay(ig,l)/(r*pt(ig,l))) + enddo + enddo + endif + + enddo ! end of nq + +end subroutine mugas_prof Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3948) +++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3949) @@ -407,4 +407,5 @@ ! --------------------------------- real gzlat(ngrid,nlayer) ! Altitude-Latitude-dependent gravity (this should be stored elsewhere...). + real pdqmuchem(ngrid,nlayer,nq) ! Condensable gases through muphi tendency (kg/kg_of_air/s). real pdqmufi(ngrid,nlayer,nq) ! Microphysical tendency (X/kg_of_air/s). real pdqmufi_prod(ngrid,nlayer,nq) ! Aerosols production tendency (kg/kg_of_air/s). @@ -1519,7 +1520,7 @@ if (tracer) then -! --------------------------------------- -! Methane ice condensation in the atmosphere -! ---------------------------------------- + !------------------------------------------------- + ! Methane (CH4) ice condensation in the atmosphere + !------------------------------------------------- rice_ch4(:,:)=0 ! initialization needed for callsedim zdqch4cloud(:,:,:)=0. @@ -1553,7 +1554,7 @@ end if -! --------------------------------------- -! CO ice condensation in the atmosphere -! ---------------------------------------- + !-------------------------------------- + ! CO ice condensation in the atmosphere + !-------------------------------------- zdqcocloud(:,:,:)=0. IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN @@ -1587,4 +1588,13 @@ rice_co(:,:)=0 ! initialization needed for callsedim END IF ! of IF (carbox) + + !-------------------------------------------------- + ! Condensable gases through muphi in the atmosphere + !-------------------------------------------------- + pdqmuchem(:,:,:) = 0. + IF (callmufi) THEN + call mugas_prof(ngrid,nlayer,nq,zzlay,zzlev,pplay,pt,pdqmuchem) + ENDIF + pdq(:,:,:) = pdq(:,:,:) + pdqmuchem(:,:,:) ! ---------------------------------------- Index: trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90 (revision 3948) +++ trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90 (revision 3949) @@ -57,6 +57,12 @@ integer,save :: igcm_n2 integer,save :: igcm_ar - integer,save :: igcm_ch4_gas ! methane gas + integer,save :: igcm_ch4_gas ! CH4 gas + integer,save :: igcm_C2H2_mugas ! C2H2 gas (condensable in microphysics) + integer,save :: igcm_C2H6_mugas ! C2H6 gas (condensable in microphysics) + integer,save :: igcm_C4H2_mugas ! C4H2 gas (condensable in microphysics) + integer,save :: igcm_C6H6_mugas ! C6H6 gas (condensable in microphysics) + integer,save :: igcm_HCN_mugas ! HCN gas (condensable in microphysics) !$OMP THREADPRIVATE(igcm_co_gas,igcm_n2,igcm_ar,igcm_ch4_gas) +!$OMP THREADPRIVATE(igcm_C2H2_mugas,igcm_C2H6_mugas,igcm_C4H2_mugas,igcm_C6H6_mugas,igcm_HCN_mugas) ! Other tracers integer,save :: igcm_ar_n2 ! for simulations using co2 + neutral gaz