Changeset 4173 for trunk

Timestamp:

Apr 7, 2026, 4:13:35 PM (4 weeks ago)

Author:

cpetetin

Message:

TITAN PCM : evaporation and small conservation correction CP

Location:

trunk/LMDZ.TITAN/libf

Files:

• muphytitan/mm_clouds.f90 (modified) (1 diff)
• muphytitan/mm_globals.f90 (modified) (1 diff)
• muphytitan/mm_haze.F90 (modified) (2 diffs)
• muphytitan/mm_microphysic.f90 (modified) (2 diffs)
• phytitan/calchim.F90 (modified) (3 diffs)
• phytitan/comchem_h.F90 (modified) (1 diff)
• phytitan/evapCH4.F90 (modified) (2 diffs)
• phytitan/physiq_mod.F90 (modified) (7 diffs)

trunk/LMDZ.TITAN/libf/muphytitan/mm_clouds.f90

@@ -391 +391 @@
     ! compute tendencies:
     ! all of these tendencies are in X.m-3 !
-    dm0aer = (zm0a - sm0a)*mm_rhoair
-    dm3aer = (zm3a - sm3a)*mm_rhoair
+    !dm0aer = (zm0a - sm0a)*mm_rhoair
+    !dm3aer = (zm3a - sm3a)*mm_rhoair
     dm0ccn = (zm0n - sm0n)*mm_rhoair
     dm3ccn = (zm3n - sm3n)*mm_rhoair
+    dm0aer = -dm0ccn
+    dm3aer = -dm3ccn
 
     dm3iX  = (zm3iX - sm3iX)                    ! this one in X.kg-1 (temporary) !

trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90

@@ -184 +184 @@
   LOGICAL, SAVE :: mm_var_prod = .false. !! Time variation of production rate control flag.
 
-  LOGICAL, SAVE :: mm_use_effg = .true. !! Enable/Disable effective G for computations.
+  LOGICAL, SAVE :: mm_use_effg = .false. !! Enable/Disable effective G for computations.
 
   !> Enable/Disable __Fiadero__'s correction.

trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.F90

@@ -773 +773 @@
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s !! Tendency of M3 (\(m^{3}.m^{-3}\)).
     INTEGER                     :: i
-    REAL(kind=mm_wp)            :: zprod,cprod,timefact
+    REAL(kind=mm_wp)            :: zprod,cprod,timefact,prod_eff
     REAL(kind=mm_wp), PARAMETER :: sigz  = 20e3_mm_wp, &
                                    fnorm = 1._mm_wp/(dsqrt(2._mm_wp*mm_pi)*sigz), &
@@ -812 +812 @@
     ENDIF
 
+    ! Correctif guassian loses when zprod is too close to border
+    prod_eff = SUM(dm3s(:)*mm_dzlev(1:mm_nla))/(mm_tx_prod*0.75_mm_wp/mm_pi * mm_dt/mm_rhoaer)
+    dm3s(:) = dm3s(:)/prod_eff
+    dm0s(:) = dm0s(:)/prod_eff
+
   END SUBROUTINE mm_haze_production
 

trunk/LMDZ.TITAN/libf/muphytitan/mm_microphysic.f90

@@ -205 +205 @@
     REAL(kind=mm_wp), INTENT(out), OPTIONAL, DIMENSION(:,:)   :: grate     !! Growth rate (\(m^{2}.s^{-1}\)).
 
-    IF (PRESENT(aer_prec))   aer_prec   = ABS(mm_aer_prec) / dt
+    IF (PRESENT(aer_prec))   aer_prec   = -mm_aer_prec / dt
     IF (PRESENT(aer_s_w))    aer_s_w    = -mm_m3as_vsed(mm_nla:1:-1)
     IF (PRESENT(aer_f_w))    aer_f_w    = -mm_m3af_vsed(mm_nla:1:-1)
@@ -212 +212 @@
 
     IF (mm_w_clouds) THEN
-      IF (PRESENT(ccn_prec))   ccn_prec   = ABS(mm_ccn_prec) / dt
-      IF (PRESENT(ice_prec))   ice_prec   = ABS(mm_ice_prec) / dt
+      IF (PRESENT(ccn_prec))   ccn_prec   = -mm_ccn_prec / dt
+      IF (PRESENT(ice_prec))   ice_prec   = -mm_ice_prec / dt
       IF (PRESENT(ccn_w))      ccn_w      = mm_ccn_vsed(mm_nla:1:-1)
       IF (PRESENT(ccn_flux))   ccn_flux   = mm_ccn_flux(mm_nla:1:-1)

trunk/LMDZ.TITAN/libf/phytitan/calchim.F90

@@ -68 +68 @@
   USE comcstfi_mod, ONLY: g, rad, pi, r, kbol
   USE geometry_mod, ONLY: latitude
+  USE mod_phys_lmdz_para, only : is_master
 #ifndef MESOSCALE  
   USE logic_mod, ONLY: moyzon_ch
@@ -166 +167 @@
   IF (firstcall) THEN
 
-     PRINT*, 'CHIMIE, premier appel'
+     if(is_master) PRINT*, 'CHIMIE, premier appel'
 
      if (ngrid .eq. 1) then ! if 1D no dynamic mixing, we set the kedd in all column
@@ -224 +225 @@
      ! a. For GCM layers we just copy-paste ( assuming that physiq always send correct altitudes ! )
 
-     PRINT*,'Init chemistry : pressure, density, temperature ... :'
-     PRINT*,'level, press_c (mbar), nb (cm-3), temp_c (K)'
+     if(is_master) then 
+         PRINT*,'Init chemistry : pressure, density, temperature ... :'
+         PRINT*,'level, press_c (mbar), nb (cm-3), temp_c (K)'
+     endif
 
 #ifndef MESOSCALE

trunk/LMDZ.TITAN/libf/phytitan/comchem_h.F90

@@ -52 +52 @@
    REAL, DIMENSION(nkim), PARAMETER               :: cmmol = (/ &
        1.01        , 2.0158      , 13.02       , 14.03       , 14.03       , 15.03       , &
-       16.04       , 24.02       , 25.03       , 26.04       , 27.05       , 28.05       , &
-       29.06       , 30.07       , 39.06       , 41.07       , 42.08       , 43.09       , &
+       16.042      , 24.02       , 25.03       , 26.036      , 27.05       , 28.05       , &
+       29.06       , 30.068      , 39.06       , 41.07       , 42.08       , 43.09       , &
        49.05       , 51.07       , 52.08       , 40.07       , 40.07       , 44.11       , &
        50.06       , 54.09       , 58.13       , 78.1136     , 38.05       , 53.07       , &
-       77.1136     , 28.0134     , 14.01       , 26.02       , 27.04       , 28.05       , &
-       39.05       , 40.04       , 41.05       , 50.04       , 51.05       , 52.04       , &
+       77.1136     , 28.0134     , 14.01       , 26.02       , 27.028      , 28.05       , &
+       39.05       , 40.04       , 41.05       , 50.04       , 51.048      , 52.04       , &
        76.1        /)
 

trunk/LMDZ.TITAN/libf/phytitan/evapCH4.F90

@@ -77 +77 @@
 !-------------
 real, parameter :: karman = 0.4 ! Karman constant [-]
-real, parameter :: humCH4 = 0.4 ! Imposed surface humidity for CH4 [-]
+real, parameter :: humCH4 = 1. ! Imposed surface humidity for CH4 [-]
 
 real, parameter :: Flnp = 0.07  ! Fraction occupied by lakes (North Pole)
@@ -129 +129 @@
 
   ! Saturation profile of CH4 [mol.mol-1] (Fray and Schmidt 2009) :
-  Tlake = tsurf(ig) - 7 ! Lakes are 2-7K less than surface.
+  Tlake = tsurf(ig) ! Lakes are 2-7K less than surface.
   qsatCH4 = (1.0e5 / pplev(ig,1)) * exp(1.051e1 - 1.110e3/Tlake - 4.341e3/Tlake**2 + 1.035e5/Tlake**3 - 7.910e5/Tlake**4)
   ! CH4 : 0.70 * qsat because of dissolution in N2
   qsatCH4 = 0.70 * qsatCH4
 
+  ! Surface humidity :
+  qsatCH4 = humCH4 * qsatCH4
+
   ! Flux at the surface [kg.m-2.s-1] :
   flux = rhoair * Cd * ws
   
-  ! Surface humidity :
-  qsatCH4 = humCH4 * qsatCH4
-
   ! <North Pole>
   if (REAL(latitude_deg(ig)) .ge. 70.0) then

trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

@@ -360 +360 @@
       ! Temporary fraction of CH4 [mol/mol]
       real, dimension(ngrid,nlayer) :: tpq_CH4
+      ! Temporary tracers tendencie [(X)/kg/s]
+      real, dimension(ngrid,nlayer,nq) :: tpdq
 
       real :: i2e(ngrid,nlayer) ! int 2 ext factor (X.kg-1 -> X.m-3)
@@ -683 +685 @@
       taux(:)        = 0.D0
       tauy(:)        = 0.D0
-      ychimbar(:,:,:)  = 0.D0
+      ychimbar(:,:,:)   = 0.D0
+      zdqmufibar(:,:,:) = 0.D0
 
       zday=pday+ptime ! Compute time, in sols (and fraction thereof).
@@ -1146 +1149 @@
             !        - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if
             !          you want to compute optics from radii.
+            tpdq(:,:,:) = pdq(:,:,:)
+            i2e(:,:) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer)
             WHERE ((pq(:,:,1)+pdq(:,:,1)*ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)*ptimestep < 0.D0))
                pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep 
@@ -1154 +1159 @@
                pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
             ENDWHERE
+            DO ig=1,ngrid    ! Add the small production to the aer loss for diagnostic balance CP
+                DO l=1, nlayer ! be carfull the rho of aerosol hardcoded 
+                        mmd_aer_prec(ig) = mmd_aer_prec(ig) - ((pdq(ig,l,2) - tpdq(ig,l,2)) + (pdq(ig,l,4) - tpdq(ig,l,4)))*i2e(g,l)*600.
+                ENDDO
+            ENDDO
             IF (callclouds) THEN
                WHERE ((pq(:,:,5)+pdq(:,:,5)*ptimestep < 0.D0) .OR. (pq(:,:,6)+pdq(:,:,6)*ptimestep < 0.D0))
@@ -1159 +1169 @@
                   pdq(:,:,6) = (epsilon(1.0)-1.D0)*pq(:,:,6)/ptimestep
                ENDWHERE
+               DO ig=1,ngrid    ! Add the small production to the ccn loss for diagnostic balance CP
+                  DO l=1, nlayer ! be carfull the rho of ccn is hardcorded
+                     mmd_ccn_prec(ig) = mmd_ccn_prec(ig) - (pdq(ig,l,6) - tpdq(ig,l,6))*i2e(ig,l)*600.
+                  ENDDO
+               ENDDO
                DO iq = 1, size(ices_indx)
                   ! For ices :
@@ -1220 +1235 @@
             ! processes ( convadj ... ) we miss them in any case as we work in zonally/diurnal 
             ! mean -> no fine diurnal/short time evolution, only seasonal evolution only.
-            if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
+            if (moyzon_ch) then
               do iq = 1,nkim
                 ychimbar(:,:,iq) =  zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
-                  if ( callclouds ) then
+                  if ( callclouds .and. mod(icount-1,ichim).eq.0 ) then
                     ychimbar(:,:,iq) =  ychimbar(:,:,iq) + ( zdqmufibar(:,:,iq+nmicro)*ptimestep / rat_mmol(iq+nmicro) )
                   endif
@@ -1319 +1334 @@
              
                ! Evaporation of lakes :
-               if (moyzon_ch) then
-                 tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
-               else
-                 tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep    ! + dycchi because ychim not yet updated [mol/mol]
-               endif
+               
+               tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep    ! + dycchi because ychim not yet updated [mol/mol]
+               
                call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
                             pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)