Index: trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3931) +++ trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3932) @@ -8,6 +8,6 @@ logical,save :: callconduct,callmolvis,callmoldiff !$OMP THREADPRIVATE(callconduct,callmolvis,callmoldiff) - logical,save :: season,diurnal,lwrite -!$OMP THREADPRIVATE(season,diurnal,lwrite) + logical,save :: season,diurnal,lwrite,evol1d +!$OMP THREADPRIVATE(season,diurnal,lwrite,evol1d) logical,save :: callgasvis,continuum,graybody !$OMP THREADPRIVATE(callgasvis,continuum,graybody) Index: trunk/LMDZ.PLUTO/libf/phypluto/evolch4.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/evolch4.F90 (revision 3932) +++ trunk/LMDZ.PLUTO/libf/phypluto/evolch4.F90 (revision 3932) @@ -0,0 +1,89 @@ + subroutine evolch4(ngrid,nlayer,pls,pplev,pdpsrf,zqch4evol) + use datafile_mod + use comcstfi_mod, only: pi + use mod_phys_lmdz_para, only : is_master + use tracer_h, only: igcm_ch4_gas, igcm_n2, mmol + + implicit none + +!================================================================== +! Purpose +! ------- +! Get tracer fields according to the solar longitude (for 1D purpose for now) +! +! Inputs +! ------ +! pls solar longitude +! pplev pressure +! pdpsrf tendency on pressure +! +! Outputs +! ------- +! zdqch4evol tendency on zqch4 + +! Both +! ---- +! +! Authors +! ------- +! Tanguy Bertrand +!================================================================== + +!----------------------------------------------------------------------- +! Arguments + + INTEGER ngrid, nlayer + REAL,INTENT(IN) :: pls + REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure levels/ + REAL,INTENT(IN) :: pdpsrf(ngrid) + REAL,INTENT(OUT) :: zqch4evol(nlayer) ! Final tendancy +!----------------------------------------------------------------------- +! Local variables + LOGICAL,SAVE :: firstcall=.true. +!$OMP THREADPRIVATE(firstcall) + + integer Nfine + parameter(Nfine=143) + integer ifine + character(len=100) :: file_path + real,save,allocatable :: lssdat(:) + real,save,allocatable :: vmrdat(:) + real :: vmrsrf + + +!---------------- INPUT ------------------------------------------------ + + IF (firstcall) then + firstcall=.false. + +!$OMP MASTER + file_path=trim(datadir)//'/cycle_vmrch4.txt' + if (is_master) print*,file_path + + open(222,file=file_path,form='formatted') + + if(.not.allocated(lssdat)) then + allocate(lssdat(Nfine)) + endif + if(.not.allocated(vmrdat)) then + allocate(vmrdat(Nfine)) + endif + + do ifine=1,Nfine + read(222,*) lssdat(ifine), vmrdat(ifine) + enddo + close(222) + +!$OMP END MASTER +!$OMP BARRIER + ENDIF + + CALL interp_line(lssdat,vmrdat,Nfine,pls*180./pi,vmrsrf,1) + + ! Reconstruct CH4 profile in percent + zqch4evol(:)=min(vmrsrf*(pplev(1,1:nlayer)/pplev(1,1))**((16.-28.)/28.),50.) + ! Convert in MMR + zqch4evol(:)=zqch4evol(:)/100.*mmol(igcm_n2)/mmol(igcm_ch4_gas) + + end subroutine evolch4 + Index: trunk/LMDZ.PLUTO/libf/phypluto/evolps.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/evolps.F90 (revision 3932) +++ trunk/LMDZ.PLUTO/libf/phypluto/evolps.F90 (revision 3932) @@ -0,0 +1,78 @@ + subroutine evolps(pls,pps) + use datafile_mod + use comcstfi_mod, only: pi + use mod_phys_lmdz_para, only : is_master + + implicit none + +!================================================================== +! Purpose +! ------- +! Get tracer fields according to the solar longitude (for 1D purpose for now) +! +! Inputs +! ------ +! pls solar longitude +! +! Outputs +! ------- +! pps surface pressure + +! Both +! ---- +! +! Authors +! ------- +! Tanguy Bertrand +!================================================================== + +!----------------------------------------------------------------------- +! Arguments + + REAL,INTENT(IN) :: pls + REAL,INTENT(OUT) :: pps +!----------------------------------------------------------------------- +! Local variables + LOGICAL,SAVE :: firstcall=.true. +!$OMP THREADPRIVATE(firstcall) + + integer Nfine + parameter(Nfine=143) + integer ifine + character(len=100) :: file_path + real,save,allocatable :: lssdat(:) + real,save,allocatable :: psdat(:) + + +!---------------- INPUT ------------------------------------------------ + + IF (firstcall) then + firstcall=.false. + +!$OMP MASTER + file_path=trim(datadir)//'/cycle_ps.txt' + if (is_master) print*,file_path + + open(222,file=file_path,form='formatted') + + if(.not.allocated(lssdat)) then + allocate(lssdat(Nfine)) + endif + if(.not.allocated(psdat)) then + allocate(psdat(Nfine)) + endif + + do ifine=1,Nfine + read(222,*) lssdat(ifine), psdat(ifine) + enddo + close(222) + +!$OMP END MASTER +!$OMP BARRIER + ENDIF + + CALL interp_line(lssdat,psdat,Nfine,pls*180./pi,pps,1) + !if (is_master) write(*,*) 'pps=',pps + + end subroutine evolps + Index: trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90 (revision 3931) +++ trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90 (revision 3932) @@ -9,6 +9,7 @@ use geometry_mod, only: longitude, latitude ! in radians use callkeys_mod, only: hazeconservch4, haze_ch4proffix, diurnal, tcon_ch4, k_ch4, & - ncratio_ch4,triton + ncratio_ch4,triton,global1d use datafile_mod, only: datadir + use write_output_mod, only: write_output implicit none @@ -23,10 +24,21 @@ ! ngrid Number of vertical columns ! nlayer Number of layers +! nq Number of tracers ! pplay(ngrid,nlayer) Pressure layers ! pplev(ngrid,nlayer+1) Pressure levels -! +! zzlay(ngrid,nlayer+1) mid layer altitude +! ptimestep Physical timestep +! zday, mu0 day clock, cos(incident flux angle) +! pq, pdq tracer MMR and tendency +! pdist_sol, declin Sun-planet distance, subsolar latitude +! pfluxuv Lyman alpha flux at Earth + ! Outputs ! ------- -! +! zdqhaze Tendency on tracers MMR haze +! zdqphot_prec Photolysis rate for haze precursors +! zdqphot_ch4 Photolysis rate for CH4 +! zdqconv_prec Tendency for haze formation and precursor destruction + ! Both ! ---- @@ -43,7 +55,4 @@ INTEGER ngrid, nlayer, nq -! REAL lati(ngrid) -! REAL long(ngrid) -! REAL declin LOGICAL firstcall SAVE firstcall @@ -55,11 +64,13 @@ REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) -! REAL,INTENT(IN) :: mmol(nq) REAL,INTENT(IN) :: pdist_sol ! distance SUN-pluto in AU - REAL,INTENT(IN) :: pfluxuv ! Lyman alpha flux at specific Ls (ph/cm2/s) + REAL,INTENT(IN) :: pfluxuv ! Lyman alpha flux at Earth at specific Ls (ph/cm2/s) REAL,INTENT(IN) :: mu0(ngrid) ! cosinus of solar incident flux - REAL,INTENT(IN) :: declin ! distance SUN-pluto in AU + REAL,INTENT(IN) :: declin ! Subsolar latitude REAL,INTENT(OUT) :: zdqhaze(ngrid,nlayer,nq) ! Final tendancy -! REAL,INTENT(OUT) :: zdqhaze_col(ngrid) ! Final tendancy haze prod + REAL,INTENT(OUT) :: zdqphot_ch4(ngrid,nlayer) ! tendancy due to photolysis of ch4 + REAL,INTENT(OUT) :: zdqphot_prec(ngrid,nlayer) ! tendancy due to photolysis of ch4 + REAL,INTENT(OUT) :: zdqconv_prec(ngrid,nlayer) ! tendancy due to conversion of + ! prec_haze into haze !----------------------------------------------------------------------- ! Local variables @@ -79,8 +90,4 @@ REAL gradflux(nlayer) ! gradient flux Lyman alpha ph.m-2.s-1 - REAL zdqphot_ch4(ngrid,nlayer) ! tendancy due to photolysis of ch4 - REAL zdqphot_prec(ngrid,nlayer) ! tendancy due to photolysis of ch4 - REAL zdqconv_prec(ngrid,nlayer) ! tendancy due to conversion of - ! prec_haze into haze REAL kch4 ! fraction of Carbon from ch4 directly dissociated ! by prec_haze @@ -104,33 +111,62 @@ REAL dist REAL long2 -!----------------------------------------------------------------------- - -!---------------- INPUT ------------------------------------------------ + + ! Diagnostics (temporary) + REAL fluxlym_sol_tot(ngrid) + REAL fluxlym_ipm_tot(ngrid) + +!----------------------------------------------------------------------- +! Input parameters avogadro = 6.022141e23 - sigch4 = 1.85e-17 ! aborption cross section ch4 cm-2 mol-1 -!! Initial Solar flux Lyman alpha on Earth (1AU) - flym_sol_earth=pfluxuv*1.e15 ! ph.m-2.s-1 -> 5e+11 ph.cm-2.s-1 -!! Initial Solar flux Lyman alpha on Pluto + sigch4 = 1.85e-17 ! aborption cross section of ch4 in cm-2 mol-1 + ! Initial Solar flux Lyman alpha at Earth (1AU) + flym_sol_earth=pfluxuv*1.e15 ! ph.m-2.s-1 -> about 5e+11 ph.cm-2.s-1 + ! Initial Solar flux Lyman alpha on Pluto flym_sol_pluto=flym_sol_earth/(pdist_sol*pdist_sol) ! ph.m-2.s-1 -! option decrease by 12% the initial IPM flux to account for Interplanetary H absorption: -! Fig 3 Gladstone et al 2014 : Lyalpha at Pluto + ! option decrease by 12% the initial Solar Lyman alpha flux to account for Interplanetary H absorption: + ! Cf Fig 3 Gladstone et al 2014 : Lyalpha at Pluto flym_sol_pluto=flym_sol_pluto*0.878 -!!!! Top of atm. Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 -! fit Fig. 4 in Randall Gladstone et al. (2015) -! -- New version : Integration over semi sphere of Gladstone data. -! Fit of results - - IF (ngrid.eq.1) THEN - flym_ipm(1)=75.e10 - mu_ipm(1) = 0.5 !max(mu0(1), 0.5) - mu_sol(1)=0.25 + ! Time constant of conversion in aerosol [second] + ! To be explored: 1.E5 - 1.E9 + tcon= tcon_ch4 + ! Parameter of conversion precurseur to haze + kch4=k_ch4 + ncratio=ncratio_ch4 ! boost for haze considering nitrogen contribution + ! ratio n/c : =0.25 if there is 1N for 3C + IF (firstcall) then + write(*,*) 'hazecloud: haze parameters:' + write(*,*) 'tcon, kch4, ncratio = ' , tcon, kch4, ncratio + firstcall=.false. + ENDIF + ! note: mu0 = cos(solar zenith angle) + ! max(mu0) at lat = declin + +!----------------------------------------------------------------------- +! Total IPM Lyman alpha flux at top of atmosphere + + ! Top of atm. Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 + ! fit Fig. 4 in Randall Gladstone et al. (2015) + ! Integration over semi sphere of Gladstone data. + ! 145 R averaged over the sky + ! 72.5 R in average over the visible hemisphere + IF (ngrid.eq.1.and..not.global1d) THEN + mu_sol(1) = mu0(1) + mu_ipm(1) = max(mu_sol(1), 0.5) + flym_ipm(1)=mu0(1)*72.5e10 + ELSE IF (ngrid.eq.1.and.global1d) THEN + mu_sol(1) = mu0(1) ! Full visible disk + mu_ipm(1) = mu0(1)/2. ! Full sphere (day+night) + flym_ipm(1)=mu_ipm(1)*145.e10 + ! Proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) ? + ! In theory, IPM flux remains relatively constant further than 15 AU + flym_ipm(1)=flym_ipm(1)*pfluxuv/5. ELSE IF (.not.diurnal) THEN - flym_ipm(:)= mu0(:)*75.e10 mu_sol(:) = mu0(:) mu_ipm(:) = max(mu_sol(:), 0.5) - ELSE -! 1) get longitude/latitude (rad) of anti-subsolar point (max de mu0 - 180) + flym_ipm(:)= mu0(:)*72.5e10 + ELSE ! case with full fit to Gladstone et al. results +! 1) get longitude/latitude (in radian) of anti-subsolar point (max de mu0 - 180) longit=-(zday-int(zday))*2.*pi IF (longit.LE.-pi) THEN @@ -143,5 +179,4 @@ valmin_dl=74.5e10 ! daylight minimum value puis=3.5 - ! 3) Loop for each location and fit DO ig=1,ngrid @@ -149,5 +184,6 @@ mu_sol(ig) = mu0(ig) mu_ipm(ig) = max(mu_sol(ig), 0.5) - IF (mu0(ig).LT.1.e-4) THEN + IF (mu0(ig).LT.1.e-4) THEN ! Daytime + ! Distance to subsolar point dist=acos(sin(latitude(ig))*sin(latit)+cos(latitude(ig))* & cos(latit)*cos(longit-longitude(ig)))*180/pi @@ -155,14 +191,16 @@ flym_ipm(ig)=(valmin_dl-valmin)/(90.**puis)*(dist)**puis & +valmin - ELSE + ELSE ! Nightime flym_ipm(ig)= mu0(ig)*(valmax-valmin_dl)+valmin_dl ENDIF - ! proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) - flym_ipm(ig)=flym_ipm(ig)*pfluxuv/5. + ! Proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) ? + ! In theory, IPM flux remains relatively constant further than 15 AU + ! flym_ipm(ig)=flym_ipm(ig)*pfluxuv/5. ENDDO ! nlayer - ENDIF ! ngrid -!--- + +!----------------------------------------------------------------------- +! Methane profile ! If fixed profile of CH4 gas @@ -179,26 +217,18 @@ close(115) endif + ! Initialization: vmrch4(:,:)=0. -!! Time constant of conversion in aerosol [second] -!! To be explore: 1.E3; 1.E5; 1.E7; 1.E9 - tcon= tcon_ch4 -!! Parameter of conversion precurseur to haze - kch4=k_ch4 - ncratio=ncratio_ch4 ! boost for haze considering nitrogen contribution - ! ratio n/c : =0.25 if there is 1N for 3C - IF (firstcall) then - write(*,*) 'hazecloud: haze parameters:' - write(*,*) 'tcon, kch4, ncratio = ' , tcon, kch4, ncratio - firstcall=.false. - ENDIF -! note: mu0 = cos(solar zenith angle) -! max(mu0) = declin - -!!---------------------------------------------------------- -!!---------------------------------------------------------- + ! Diagnostics (temporary) + fluxlym_sol_tot(:)=flym_sol_pluto*mu_sol(:) + fluxlym_ipm_tot(:)=flym_ipm(:) + call write_output("fluxlym_sol","solar lyman alpha flux","",fluxlym_sol_tot) + call write_output("fluxlym_ipm","IPM lyman alpha flux","",fluxlym_ipm_tot) +!----------------------------------------------------------------------- +! Main Loop DO ig=1,ngrid + !! Initializations zq_ch4(ig,:)=0. zq_prec(ig,:)=0. @@ -210,8 +240,9 @@ gradflux(:)=0. fluxlym_sol(1:nlayer)=0. + fluxlym_sol(nlayer+1)=flym_sol_pluto*mu_sol(ig) ! fluxlym_ipm(nlayer+1)= flym_ipm(ig) - !! Interpolate on the model vertical grid + !! Interpolate CH4 MMR on the model vertical grid if (haze_ch4proffix) then CALL interp_line(levdat,vmrdat,Nfine, & @@ -220,4 +251,5 @@ DO l=nlayer,1,-1 + !! Actualisation tracer ch4 and prec_haze if (haze_ch4proffix) then @@ -229,15 +261,17 @@ endif + zq_prec(ig,l)=pq(ig,l,igcm_prec_haze)+ & + pdq(ig,l,igcm_prec_haze)*ptimestep + + !! Sanity check if (zq_ch4(ig,l).lt.0.) then zq_ch4(ig,l)=0. endif - zq_prec(ig,l)=pq(ig,l,igcm_prec_haze)+ & - pdq(ig,l,igcm_prec_haze)*ptimestep - - !! Calculation optical depth ch4 in Lyman alpha for each layer l + !! Calculation of CH4 optical depth in Lyman alpha for each layer l tauch4(l)=sigch4*1.e-4*avogadro* & (zq_ch4(ig,l)/(mmol(igcm_ch4_gas)*1.e-3))* & (pplev(ig,l)-pplev(ig,l+1))/g + !! Calculation of Flux in each layer l if (mu_sol(ig).gt.1.e-6) then @@ -249,4 +283,5 @@ gradflux(l)=fluxlym_sol(l+1)-fluxlym_sol(l) & + fluxlym_ipm(l+1)-fluxlym_ipm(l) + !! tendancy on ch4 !! considering 1 photon destroys 1 ch4 by photolysis @@ -254,8 +289,8 @@ gradflux(l)*g/(pplev(ig,l)-pplev(ig,l+1)) - !! tendency of prec created by photolysis of ch4 + !! tendency of precursors created by photolysis of ch4 zdqphot_prec(ig,l)=-zdqphot_ch4(ig,l) - !! update precurseur zq_prec + !! update of precursors mass: zq_prec zq_prec(ig,l)=zq_prec(ig,l)+ & zdqphot_prec(ig,l)*ptimestep @@ -269,9 +304,7 @@ ENDDO ! nlayer - !! Final tendancy for prec_haze and haze - + !! Final tendancy for prec_haze, haze and CH4 DO iq=1,nq tname=noms(iq) - !print*, 'TB17 tname=',tname,tname(1:4) if (tname(1:4).eq."haze") then zdqhaze(ig,:,iq) = -zdqconv_prec(ig,:)* & @@ -280,5 +313,4 @@ zdqhaze(ig,:,igcm_prec_haze)= zdqphot_prec(ig,:) + & zdqconv_prec(ig,:) - else if (noms(iq).eq."ch4_gas".and.(.not.hazeconservch4).and.(.not.haze_ch4proffix)) then zdqhaze(ig,:,igcm_ch4_gas)= zdqphot_ch4(ig,:) @@ -288,14 +320,4 @@ ENDDO ! ngrid - - !! tendency kg/m2/s for haze column: -! zdqhaze_col(:)=0. -! DO ig=1,ngrid -! DO l=1,nlayer -! zdqhaze_col(ig)=zdqhaze_col(ig)+zdqhaze(ig,l,igcm_haze)* & -! (pplev(ig,l)-pplev(ig,l+1))/g -! ENDDO -! ENDDO - end subroutine hazecloud Index: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3931) +++ trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3932) @@ -213,4 +213,8 @@ if (is_master) write(*,*) trim(rname)//": season = ",season + evol1d=.false. ! default value + call getin_p("evol1d",evol1d) + if (is_master) write(*,*) trim(rname)//": evol1d = ",evol1d + if (is_master) then write(*,*) trim(rname)//": Fast mode (nogcm) ?" @@ -1351,4 +1355,7 @@ call getin_p("kmixmin",kmixmin) if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin + kmix_proffix=.false. ! default value + call getin_p("kmix_proffix",kmix_proffix) + if (is_master) write(*,*)trim(rname)//": kmix_proffix = ",kmix_proffix if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1' Index: trunk/LMDZ.PLUTO/libf/phypluto/lymalpha.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/lymalpha.F90 (revision 3931) +++ trunk/LMDZ.PLUTO/libf/phypluto/lymalpha.F90 (revision 3932) @@ -17,4 +17,5 @@ ! Outputs ! ------- +! pflux Lyman alpha flux ! ! Both @@ -33,6 +34,4 @@ !----------------------------------------------------------------------- ! Local variables - REAL :: vectls - REAL :: vectflux LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) @@ -40,5 +39,5 @@ !!read lyman alpha flux integer Nfine - parameter(Nfine=13281) + parameter(Nfine=9047) integer ifine character(len=100) :: file_path Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3931) +++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3932) @@ -66,5 +66,5 @@ tracer, UseTurbDiff, & global1d, szangle, & - callmufi + callmufi, evol1d use check_fields_mod, only: check_physics_fields use conc_mod, only: rnew, cpnew, ini_conc_mod @@ -362,4 +362,6 @@ real zdqmoldiff(ngrid,nlayer,nq) + REAL zqch4evol(nlayer) + ! Haze relatated tendancies REAL zdqhaze(ngrid,nlayer,nq) @@ -379,4 +381,5 @@ REAL flym_ipm(ngrid) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 REAL zfluxuv ! Lyman alpha flux at 1AU + REAL pss ! evolving surface pressure in 1D REAL array_zero1(ngrid) @@ -736,5 +739,5 @@ end if -! Get Lyman alpha flux at specific Ls + ! Get Lyman alpha flux at specific Ls if (callmufi.or.haze) then call lymalpha(zls,zfluxuv) @@ -1438,4 +1441,22 @@ ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond') endif + endif + + ! Apply empiric surface pressure change in 1D when + ! evol1D is active (seasonal evolution of pressure) + if (ngrid.eq.1.and.season.and.evol1d) then + ! Get empiric value of surface pressure depending on Ls + call evolps(zls,pss) + ! Update surface pressure tendancy + pdpsrf(1) = pdpsrf(1)+ (pss-(pplev(1,1)+ & + pdpsrf(1)*ptimestep))/ptimestep + if (methane) then + ! Get empiric value of CH4 MMR depending on Ls + call evolch4(ngrid,nlayer,zls,pplev,pdpsrf,zqch4evol) + ! Update tracer tendency + pdq(1,:,igcm_ch4_gas)=pdq(1,:,igcm_ch4_gas)+ & + (zqch4evol(:)-(pq(1,:,igcm_ch4_gas)+ & + pdq(1,:,igcm_ch4_gas)*ptimestep))/ptimestep + endif endif