Changeset 3932 for trunk

Timestamp:

Oct 23, 2025, 11:51:31 AM (3 months ago)

Author:

tbertrand

Message:

PLUTO PCM: Allow 1D model to be run with empirical seasonal cycles of N2 and CH4 (change in pressure and CH4 abundance with time)
TB

Location:

trunk/LMDZ.PLUTO/libf/phypluto

Files:

• callkeys_mod.F90 (modified) (1 diff)
• evolch4.F90 (added)
• evolps.F90 (added)
• hazecloud.F90 (modified) (18 diffs)
• inifis_mod.F90 (modified) (2 diffs)
• lymalpha.F90 (modified) (3 diffs)
• physiq_mod.F90 (modified) (5 diffs)

trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90

@@ -8 +8 @@
       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)

trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90

@@ -9 +9 @@
       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 +24 @@
 !     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 +55 @@
 
       INTEGER ngrid, nlayer, nq
-!      REAL lati(ngrid)
-!      REAL long(ngrid)
-!      REAL declin
       LOGICAL firstcall
       SAVE firstcall
@@ -55 +64 @@
       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 +90 @@
 
       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 +111 @@
       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 +179 @@
         valmin_dl=74.5e10  ! daylight minimum value
         puis=3.5
-
 !     3) Loop for each location and fit
         DO ig=1,ngrid
@@ -149 +184 @@
           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 +191 @@
            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 +217 @@
         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 +240 @@
         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 +251 @@
 
         DO l=nlayer,1,-1
+
           !! Actualisation tracer ch4 and prec_haze
           if (haze_ch4proffix) then
@@ -229 +261 @@
           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 +283 @@
           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 +289 @@
             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 +304 @@
         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 +313 @@
               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 +320 @@
       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
 

trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90

@@ -213 +213 @@
      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 +1355 @@
      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'

trunk/LMDZ.PLUTO/libf/phypluto/lymalpha.F90

@@ -17 +17 @@
 !     Outputs
 !     -------
+!     pflux               Lyman alpha flux
 !
 !     Both
@@ -33 +34 @@
 !-----------------------------------------------------------------------
 !     Local variables
-      REAL :: vectls
-      REAL :: vectflux
       LOGICAL,SAVE :: firstcall=.true.
 !$OMP THREADPRIVATE(firstcall)
@@ -40 +39 @@
       !!read lyman alpha flux
       integer Nfine
-      parameter(Nfine=13281)
+      parameter(Nfine=9047)
       integer ifine
       character(len=100) :: file_path

trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

@@ -66 +66 @@
                               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 +362 @@
       real zdqmoldiff(ngrid,nlayer,nq)
 
+      REAL zqch4evol(nlayer)
+
       ! Haze relatated tendancies
       REAL zdqhaze(ngrid,nlayer,nq)
@@ -379 +381 @@
       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 +739 @@
       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 +1441 @@
                    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