Changeset 4146 for trunk/LMDZ.GENERIC/libf/phygeneric/inifis_mod.F90

Timestamp:

Mar 19, 2026, 2:35:46 PM (11 days ago)

Author:

gmilcareck

Message:

Thermodynamics update on LMDZ.GENERIC

File:

• trunk/LMDZ.GENERIC/libf/phygeneric/inifis_mod.F90 (modified) (5 diffs)

trunk/LMDZ.GENERIC/libf/phygeneric/inifis_mod.F90

@@ -19 +19 @@
   use time_phylmdz_mod, only: diagfi_output_rate, slow_diagfi, &
                               init_time, daysec, dtphys
-  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
-                          mugaz, pi, avocado
+  use comcstfi_mod, only: rad, cppd_ref, cppv_ref, cppc_ref, g, rd_ref, rcp_ref, &
+                          mugaz_ref, pi, avocado
   use planete_mod, only: nres
   use planetwide_mod, only: planetwide_sumval
@@ -90 +90 @@
   ! initialize constants in comcstfi_mod
   rad=prad
-  cpp=pcpp
+  cppd_ref=pcpp
   g=pg
-  r=pr
-  rcp=r/cpp
-  mugaz=8.314463*1000./pr ! dummy init
+  rd_ref=pr
+  rcp_ref=rd_ref/cppd_ref
+  mugaz_ref=8.314463*1000./pr ! dummy init
   pi=2.*asin(1.)
   avocado = 6.022141e23   ! added by RW
@@ -1103 +1103 @@
      if (is_master) write(*,*)trim(rname)//": generic_condensation = ",generic_condensation
      
+     if (generic_condensation) then
+                if (is_master) write(*,*)trim(rname)//": value for metallicity? "
+                metallicity=0.0 ! default value
+                call getin_p("metallicity",metallicity)
+                if (is_master) write(*,*)trim(rname)//": metallicity = ",metallicity
+
+                if (is_master) write(*,*)trim(rname)//": Deep generic tracer vapor mixing ratio ? (no effect if negative) "
+                qvap_deep=-1. ! default value
+                call getin_p("qvap_deep",qvap_deep)
+                if (is_master) write(*,*)trim(rname)//": qvap_deep = ",qvap_deep   
+
+                if (is_master) write(*,*)trim(rname)//":top generic tracer vapor mixing ratio ? (no effect if negative) "
+                qvap_top=-1. ! default value
+                call getin_p("qvap_top",qvap_top)
+                if (is_master) write(*,*)trim(rname)//": qvap_top = ",qvap_top 
+
+                if (is_master) write(*,*)trim(rname)//": align_strato_cold_trap ? "
+                align_strato_cold_trap=.false. ! default value
+                call getin_p("align_strato_cold_trap",align_strato_cold_trap)
+                if (is_master) write(*,*)trim(rname)//": align_strato_cold_trap = ",align_strato_cold_trap
+     endif
+     
      if (is_master) write(*,*)trim(rname)//": Generic rain of tracers ?"
      generic_rain=.false. !default value 
      call getin_p("generic_rain",generic_rain)
      if (is_master) write(*,*)trim(rname)//": generic_rain = ",generic_rain
+     
+     if (generic_rain) then
+            if (is_master) write(*,*)trim(rname)//": re-evaporate precipitations?"
+            evap_prec_generic=.true. ! default value
+            call getin_p("evap_prec_generic",evap_prec_generic)
+            if (is_master) write(*,*)trim(rname)//": evap_prec_generic = ",evap_prec_generic
+
+            if (evap_prec_generic) then
+               if (is_master) write(*,*)trim(rname)//": multiplicative constant in reevaporation"
+               evap_coeff_generic=1.   ! default value
+               call getin_p("evap_coeff_generic",evap_coeff_generic)
+               if (is_master) write(*,*)trim(rname)//": evap_coeff_generic = ",evap_coeff_generic       
+            end if
+
+            if (is_master) write(*,*)trim(rname)//": Precipitation scheme to use?"
+            precip_scheme_generic=1 ! default value
+            call getin_p("precip_scheme_generic",precip_scheme_generic)
+            if (is_master) write(*,*)trim(rname)//": precip_scheme_generic = ",precip_scheme_generic
+
+            if (precip_scheme_generic.eq.1) then
+               if (is_master) write(*,*)trim(rname)//": rainthreshold_generic in simple scheme?"
+               rainthreshold_generic=0. ! default value
+               call getin_p("rainthreshold_generic",rainthreshold_generic)
+               if (is_master) write(*,*)trim(rname)//": rainthreshold_generic = ",rainthreshold_generic
+               
+            else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then
+               
+               if (is_master) write(*,*)trim(rname)//": cloud GCS saturation level in non simple scheme?"
+               cloud_sat_generic=2.6e-4   ! default value
+               call getin_p("cloud_sat_generic",cloud_sat_generic)
+               if (is_master) write(*,*)trim(rname)//": cloud_sat_generic = ",cloud_sat_generic
+            
+               if (is_master) write(*,*)trim(rname)//": precipitation timescale in non simple scheme?"
+               precip_timescale_generic=3600.  ! default value
+               call getin_p("precip_timescale_generic",precip_timescale_generic)
+               if (is_master) write(*,*)trim(rname)//": precip_timescale_generic = ",precip_timescale_generic
+
+            else if (precip_scheme_generic.eq.4) then
+               
+               if (is_master) write(*,*)trim(rname)//": multiplicative constant in Boucher 95 precip scheme"
+               Cboucher_generic=1.   ! default value
+               call getin_p("Cboucher_generic",Cboucher_generic)
+               if (is_master) write(*,*)trim(rname)//": Cboucher_generic = ",Cboucher_generic   
+
+            endif
+     endif
 
      if (is_master) write(*,*)trim(rname)//": Include moist adjustement for GCS ?"
@@ -1207 +1275 @@
        call getin_p("nvarlayer",nvarlayer)
      endif
-
-     if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
-     if (is_master) write(*,*)'Predefined Mg from dynamics is ',mugaz,'amu'
-
-     force_cpp=.false. ! default value
-     call getin_p("force_cpp",force_cpp)
-     if (force_cpp) then
-      if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp
-      if (is_master) write(*,*)trim(rname)//": force_cpp is deprecated.",&
-      "Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
-      call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
-     endif
-
-     if (is_master) write(*,*)trim(rname)//&
-     ": where do you want your cpp/mugaz value to come from?",&
-     "=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
+     
+     if (is_master) write(*,*)trim(rname)//&
+     ": where do you want your cppd_ref/mugaz_ref value to come from?",&
+     "=> 0: dynamics (3d), 1: forced in callphys.def (1d), 2: computed from gases.def (1d)?"
      cpp_mugaz_mode = 0 ! default value
      call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
@@ -1235 +1291 @@
         call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
      endif
-
-     if (cpp_mugaz_mode == 1) then
-       mugaz = -99999.
+     
+     if (is_master) write(*,*)trim(rname)//&
+     ": CHOOSE THERMODYNAMICS MOD"
+     thermo_phy='thermo_uni_ideal' !default value 
+     call getin_p("thermo_phy",thermo_phy)
+     if (is_master) write(*,*)trim(rname)//": thermo_phy = ",trim(thermo_phy)
+     if(thermo_phy.eq.'thermo_uni_ideal') then
+       if (cpp_mugaz_mode == 0) write(*,*) "cppd_ref and mugaz_ref fixed by the dynamics"
+       if (cpp_mugaz_mode == 1) then
+         write(*,*) "cppd_ref and mugaz_ref fixed by the physics (callphys.def)"
+         cppd_ref = -99999.
+         if (is_master) write(*,*)trim(rname)//&
+           ": DRY SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
+         call getin_p("cppd_ref",cppd_ref)
+         IF (cppd_ref.eq.-99999.) THEN
+             PRINT *, "cppd_ref must be set if thermo_phy = thermo_uni_ideal and cpp_mugaz_mode == 1"
+             STOP
+         ENDIF
+         mugaz_ref = -99999.
+         if (is_master) write(*,*)trim(rname)//&
+           ": MEAN MOLECULAR MASS in g mol-1 ?"
+         call getin_p("mugaz_ref",mugaz_ref)
+         IF (mugaz_ref.eq.-99999.) THEN
+           call abort_physic(rname,"mugaz_ref must be set if cpp_mugaz_mode = 1",1)
+         ENDIF
+       endif
+       if (cpp_mugaz_mode == 2) write(*,*) "cppd_ref and mugaz_ref calculated by calc_cpp_mugaz in physics"
+       call su_gases
+       call calc_cpp_mugaz
+       if (generic_condensation) then
+         cppv_ref = -99999.
+         if (is_master) write(*,*)trim(rname)//&
+           ": VAPOUR TRACER SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
+         call getin_p("cppv_ref",cppv_ref)
+         IF (cppv_ref.eq.-99999.) THEN
+             PRINT *, "cppv_ref must be set if generic_condensation = true"
+             STOP
+         ENDIF
+       endif
+     elseif(thermo_phy.eq.'thermo_binary_ideal') then
+       write(*,*) "WARNING: thermo_binary_ideal can be used only wit the DYNAMICO CORE or in 1D"
+       if (.not.tracer) then
+          call abort_physic(rname,"Error: We need a tracer for the thermo_binary_ideal mod!",1)
+       endif 
+       if (cpp_mugaz_mode == 0) write(*,*) "cppd_ref and mugaz_ref fixed by the dynamics"
+       if (cpp_mugaz_mode == 1) then
+         write(*,*) "cppd_ref and mugaz_ref fixed by the physics (callphys.def)"
+         cppd_ref = -99999.
+         if (is_master) write(*,*)trim(rname)//&
+           ": DRY SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
+         call getin_p("cppd_ref",cppd_ref)
+         IF (cppd_ref.eq.-99999.) THEN
+             PRINT *, "cppd_ref must be set if thermo_phy = thermo_binary_ideal"
+             STOP
+         ENDIF
+         mugaz_ref = -99999.
+         if (is_master) write(*,*)trim(rname)//&
+           ": MEAN MOLECULAR MASS in g mol-1 ?"
+         call getin_p("mugaz_ref",mugaz_ref)
+         IF (mugaz_ref.eq.-99999.) THEN
+           call abort_physic(rname,"mugaz_ref must be set if cpp_mugaz_mode = 1",1)
+         ENDIF
+       endif
+       if (cpp_mugaz_mode == 2) write(*,*) "cppd_ref and mugaz_ref calculated by calc_cpp_mugaz in physics"
+       call su_gases
+       call calc_cpp_mugaz
+       cppv_ref = -99999.
        if (is_master) write(*,*)trim(rname)//&
-         ": MEAN MOLECULAR MASS in g mol-1 ?"
-       call getin_p("mugaz",mugaz)
-       IF (mugaz.eq.-99999.) THEN
-         call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
+         ": VAPOUR TRACER SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
+       call getin_p("cppv_ref",cppv_ref)
+       IF (cppv_ref.eq.-99999.) THEN
+           PRINT *, "cppv_ref must be set if thermo_phy = thermo_binary_ideal"
+           STOP
        ENDIF
-       cpp = -99999.
+       cppc_ref = -99999.
        if (is_master) write(*,*)trim(rname)//&
-         ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
-       call getin_p("cpp",cpp)
-       IF (cpp.eq.-99999.) THEN
-           call abort_physic(rname, "cpp must be set if cpp_mugaz_mode = 1", 1)
+         ": CONDENSED TRACER SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
+       call getin_p("cppc_ref",cppc_ref)
+       IF (cppc_ref.eq.-99999.) THEN
+           PRINT *, "cppc_ref must be set if thermo_phy = thermo_binary_ideal"
+           STOP
        ENDIF
-       if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
-       if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu'
-  
-     endif ! of if (cpp_mugaz_mode == 1)
-     call su_gases
-     call calc_cpp_mugaz
+     else
+       call abort_physic(rname,"thermodynamics mod not recognized",1)
+     endif
 
      if (is_master) then