Changeset 4146 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Mar 19, 2026, 2:35:46 PM (10 days ago)
Author:
gmilcareck
Message:

Thermodynamics update on LMDZ.GENERIC

Location:
trunk/LMDZ.GENERIC
Files:
1 added
1 deleted
50 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/changelog.txt

    r4137 r4146  
    22162216== 17/03/2026 == GM
    22172217The nlines were not reset to 0 when reading tracer profiles during the newstart.
     2218
     2219== 19/03/2026 == GM
     2220***THERMODYNAMICS UPDATE ON LMDZ.GENERIC ***
     2221- The specific heat capacity, the specific gas constant and the R/Cp factor
     2222  are now calculated in the new thermo_mod module.
     2223- Calculations of the potential temperature and the Exner function are now performed
     2224  in the thermodynamics subroutine. This routine is called in physiq_mod whenever
     2225  the temperature is updated after a subroutine is called.
     2226  PLEASE NOTE: for future updates to LMDZ_GENERIC, if you wish to add a new subroutine
     2227  that modifies the temperature or pressure in `physiq_mod`, you must call this
     2228  subroutine to update `theta` and `exner` (see examples on physiq_mod),
     2229  otherwise, the thermodynamics will be incorrect.
     2230- IMPORTANT: The constants cpp and mugaz are now called cppd_ref and mugaz_ref.
     2231  You must therefore update your callphys.def. The cpp_mugaz_mode option is still
     2232  valid, where mode 0 reads the cpp from the dynamics, mode 1 reads cppd_ref and
     2233  mugaz_ref from callphys.def, and mode 2 calculates cppd_ref and mugaz_ref itself
     2234  using the calc_cpp_mugaz routine.
     2235- For the moment, only the ideal single-gas option (thermo_uni_ideal) is available
     2236  (i.e. the current case). It is set by default in the model with the thermo_phy flag.
     2237  You do not need to change or add anything to your .def files.
     2238  The ideal gas mixture option will be added at a later date.
     2239- If the `generic_condensation` option is enabled, you must now add
     2240  the constant `cppv_ref`, which represents the specific heat capacity of
     2241  the pure condensable gas you have chosen.
     2242- Normally, this commit will not change your previous simulations, but,
     2243  a difference of around 1K may occur in some configurations because theta was
     2244  not well updated before for example.
     2245 === RESUME ===
     2246- If you add a new subroutine on physiq_mod that changes the temperature or pressure,
     2247  call thermodynamics after that to update potential temperature and exner function.
     2248- cpp becomes cppd_ref in callphys.def
     2249- mugaz becomes mugaz_ref in callphys.def
     2250- NEW -> add cppv_ref if generic_condensation=.true.
     2251This update changes a lot of things, then some problems can occur.
     2252Some fixes will be added as the model is tested in various configurations.
     2253

  • trunk/LMDZ.GENERIC/libf/aeronogeneric/calchim_asis.F90

    r3893 r4146  
    262262            zdtchim_output(ig) = 0.
    263263            do l = 1,nlayer
    264               dEzchim(ig,l) = zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
    265               zdtchim_output(ig) = zdtchim_output(ig) + zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
     264              dEzchim(ig,l) = zdtchim(ig,l)*cppd_ref*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
     265              zdtchim_output(ig) = zdtchim_output(ig) + zdtchim(ig,l)*cppd_ref*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
    266266            end do
    267267

  • trunk/LMDZ.GENERIC/libf/aeronogeneric/iniwrite_spec.F

    r3309 r4146  
    22
    33      use photolysis_mod, only : nw, wl, wc, wu
    4       use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
     4      use comcstfi_mod, only: rad, omeg, g, mugaz_ref, rcp_ref, pi
    55      use time_phylmdz_mod, only: daysec, dtphys
    66!      USE logic_mod, ONLY: fxyhypb,ysinus
    … …  
    7272      tab_cntrl(6)  = omeg
    7373      tab_cntrl(7)  = g
    74       tab_cntrl(8)  = mugaz
    75       tab_cntrl(9)  = rcp
     74      tab_cntrl(8)  = mugaz_ref
     75      tab_cntrl(9)  = rcp_ref
    7676      tab_cntrl(10) = daysec
    7777      tab_cntrl(11) = dtphys

  • trunk/LMDZ.GENERIC/libf/aeronogeneric/photochemistry_asis.F90

    r3893 r4146  
    1919
    2020      use callkeys_mod
    21       use comcstfi_mod, only:  r,cpp,avocado,pi
     21      use comcstfi_mod, only:  rd_ref,cppd_ref,avocado,pi
    2222      use tracer_h
    2323      use types_asis
    … …  
    203203
    204204zdtchim(:) = 0.
    205 rho(:)     = (press(:)*1e2)/(r*temp(:))
     205rho(:)     = (press(:)*1e2)/(rd_ref*temp(:))
    206206
    207207!===================================================================
    … …  
    212212
    213213  do iphot = 1,nb_phot_hv_max
    214     zdtchim(:nlayer-1) = zdtchim(:nlayer-1) + e_phot(:nlayer-1,iphot)*c(:nlayer-1,indexchim(trim(jlabel(iphot,2))))/(cpp*(rho(:nlayer-1)*1e-6)*avocado)
     214    zdtchim(:nlayer-1) = zdtchim(:nlayer-1) + e_phot(:nlayer-1,iphot)*c(:nlayer-1,indexchim(trim(jlabel(iphot,2))))/(cppd_ref*(rho(:nlayer-1)*1e-6)*avocado)
    215215  end do
    216216

  • trunk/LMDZ.GENERIC/libf/aeronogeneric/photolysis_asis.F90

    r2542 r4146  
    166166      do l = 1,lswitch-1
    167167!         colo3min    = col(l)*7.171e-10
    168          colo3min    = col(l)*1.227e-10*(g/9.81)*(mugaz/28)
     168         colo3min    = col(l)*1.227e-10*(g/9.81)*(mugaz_ref/28)
    169169         ratio_o3(l) = colo3(l)/colo3min
    170170         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo))

  • trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phygeneric/inichim_newstart.F90

    r4137 r4146  
    131131         else if (qxf(iq) /= 'None') then
    132132            ! Opening file
     133            nlines=0
    133134            fil = trim(datadir)//'/chemical_profiles/'//qxf(iq)
    134135            print*, 'chemical pofile '//trim(noms(iq))//': ', fil

  • trunk/LMDZ.GENERIC/libf/phygeneric/aerosol_opacity.F90

    r4077 r4146  
    1717                              iaero_venus3, iaero_venusUV
    1818       USE tracer_h, only: noms,rho_co2,rho_ice,rho_q,mmol
    19        use comcstfi_mod, only: g, pi, mugaz, avocado
     19       use comcstfi_mod, only: g, pi, avocado
    2020       use geometry_mod, only: latitude
    2121       use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl, &

  • trunk/LMDZ.GENERIC/libf/phygeneric/calc_cpp_mugaz.F90

    r3641 r4146  
    1919
    2020      use gases_h
    21       use comcstfi_mod, only: cpp, mugaz
     21      use comcstfi_mod, only: cppd_ref, mugaz_ref
    2222      use callkeys_mod, only: cpp_mugaz_mode
    2323      use mod_phys_lmdz_para, only : is_master
    … …  
    155155      endif
    156156
    157       if(((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
    158            (abs(1.-mugaz/mugaz_c).gt.1.e-6)).and. is_master ) then
     157      if(((abs(1.-cppd_ref/cpp_c).gt.1.e-6) .or.  &
     158           (abs(1.-mugaz_ref/mugaz_c).gt.1.e-6)).and. is_master ) then
    159159         ! Ehouarn: tolerate a small mismatch between computed/stored values
    160          print*,'--> Values do not match with the predefined one ! (',cpp,',',mugaz,')'
     160         print*,'--> Values do not match with the predefined one ! (',cppd_ref,',',mugaz_ref,')'
    161161         print*,'    Because cp varies with temperature and that some gases may not appear in gases.def,'
    162162         print*,'    a small discrepancy might be completely normal.'
    … …  
    168168      if (cpp_mugaz_mode==2) then
    169169          if (is_master) print*,'*** cpp_mugaz_mode==2, so setting cpp & mugaz to computations in calc_cpp_mugaz.'
    170           mugaz = mugaz_c
    171           cpp = cpp_c
     170          mugaz_ref = mugaz_c
     171          cppd_ref = cpp_c
    172172      else
    173173          if (is_master) print*,'*** Leaving cpp & mugaz equal to predefined values'

  • trunk/LMDZ.GENERIC/libf/phygeneric/calcul_fluxs_mod.F90

    r3100 r4146  
    1414
    1515    USE dimphy
    16     USE comcstfi_mod, ONLY: r
     16    USE comcstfi_mod, ONLY: rd_ref
    1717!    INCLUDE "YOMCST.h"
    1818!    INCLUDE "clesphys.h"
    … …  
    4444    DO i=1,knon
    4545       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
    46        buf = cdrag_m(i) * mod_wind * p1lay(i)/(r*t1lay(i))
     46       buf = cdrag_m(i) * mod_wind * p1lay(i)/(rd_ref*t1lay(i))
    4747       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
    4848       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )

  • trunk/LMDZ.GENERIC/libf/phygeneric/callkeys_mod.F90

    r4055 r4146  
    5959      logical,save :: generic_rain
    6060      logical,save :: virtual_theta_correction
    61 !$OMP THREADPRIVATE(varactive,varfixed,sedimentation,generic_condensation,generic_rain,virtual_theta_correction)
     61      character(64),save :: thermo_phy
     62!$OMP THREADPRIVATE(varactive,varfixed,sedimentation,generic_condensation,generic_rain,virtual_theta_correction,thermo_phy)
    6263      logical,save :: water ,watercond, waterrain, moistadjustment, moistadjustment_generic, moist_convection_inhibition
    6364!$OMP THREADPRIVATE(water, watercond, waterrain, moistadjustment, moistadjustment_generic, moist_convection_inhibition)
     65
     66      real, save :: metallicity !metallicity of planet
     67      real, save :: qvap_deep   ! deep mixing ratio of vapor when simulating bottom less planets
     68      real, save :: qvap_top   ! top mixing ratio of vapor when simulating bottom less planets
     69      logical, save :: align_strato_cold_trap
     70!$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, align_strato_cold_trap)
     71
     72      logical,save :: evap_prec_generic ! Does the rain evaporate ?
     73      real,save :: evap_coeff_generic ! multiplication evaporation constant. 1. gives the model of gregory et al.
     74      integer,save :: precip_scheme_generic      ! id number for precipitaion scheme
     75      real,save :: rainthreshold_generic                ! Precipitation threshold in simple scheme
     76      real,save :: cloud_sat_generic                    ! Precipitation threshold in non simple scheme
     77      real,save :: precip_timescale_generic             ! Precipitation timescale
     78      real,save :: Cboucher_generic ! Precipitation constant in Boucher 95 scheme
     79!$OMP THREADPRIVATE(evap_prec_generic,evap_coeff_generic,precip_scheme_generic,rainthreshold_generic,cloud_sat_generic,precip_timescale_generic,Cboucher_generic)
     80
     81
    6482      logical,save :: aeroco2, aeroh2o, aeroh2so4, aeroback2lay
    6583!$OMP THREADPRIVATE(aeroco2, aeroh2o, aeroh2so4, aeroback2lay)

  • trunk/LMDZ.GENERIC/libf/phygeneric/comcstfi_mod.F90

    r3663 r4146  
    55      REAL,SAVE :: rad ! radius of the planet (m)
    66      REAL,SAVE :: g ! gravity (m/s2)
    7       REAL,SAVE :: r ! reduced gas constant (r=8.314463/(mugaz/1000.0))
    8       REAL,SAVE :: cpp ! Cp of the atmosphere
    9       REAL,SAVE :: rcp ! r/cpp
    10       REAL,SAVE :: mugaz ! molar mass of the atmosphere (g/mol)
     7      REAL,SAVE :: rd_ref ! reduced dry gas constant (r=8.314463/(mugaz/1000.0))
     8      REAL,SAVE :: cppd_ref ! Total cp of the non condensable species
     9      REAL,SAVE :: cppv_ref ! cp of the tracer
     10      REAL,SAVE :: cppc_ref ! cp of the condensed tracer (liquid or solid)
     11      REAL,SAVE :: rcp_ref ! constant rd_ref/cppd_ref
     12      REAL,SAVE :: mugaz_ref ! molar mass of the dry atmosphere (g/mol)
    1113      REAL,SAVE :: omeg ! planet rotation rate (rad/s)
    12       REAL,SAVE :: avocado !  6.022141e23
    13 !$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,mugaz,omeg,avocado)
     14      REAL,SAVE :: avocado ! 6.022141e23
     15!$OMP THREADPRIVATE(pi,rad,g,rd_ref,cppd_ref,cppv_ref,cppc_ref,rcp_ref,mugaz_ref,omeg,avocado)
    1416
    1517END MODULE comcstfi_mod

  • trunk/LMDZ.GENERIC/libf/phygeneric/condensation_generic_mod.F90

    r3280 r4146  
    66    subroutine condensation_generic(ngrid,nlayer,nq,ptimestep, pplev, pplay,    &
    77                pt, pq, pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
    8         use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity
     8        use callkeys_mod, only: thermo_phy,metallicity,qvap_deep,qvap_top,align_strato_cold_trap
    99        use generic_cloud_common_h
    1010        USE tracer_h
    1111        USE mod_phys_lmdz_para, only: is_master
    1212        use generic_tracer_index_mod, only: generic_tracer_index
     13        use thermo_mod
     14        use comcstfi_mod
    1315        IMPLICIT none
    1416
    … …  
    4547        REAL, intent(out) :: rneb(ngrid,nlayer,nq)    ! fraction nuageuse
    4648
    47 !       Options :
    48         real, save :: metallicity !metallicity of planet
    49         REAL, SAVE :: qvap_deep   ! deep mixing ratio of vapor when simulating bottom less planets
    50         REAL, SAVE :: qvap_top   ! top mixing ratio of vapor when simulating bottom less planets
    51         logical, save :: align_strato_cold_trap
    52 !$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, align_strato_cold_trap)
    53 
    5449!       Local variables
    5550
    … …  
    6661        DOUBLE PRECISION zq(ngrid)
    6762        DOUBLE PRECISION zcond(ngrid),zcond_iter
    68         DOUBLE PRECISION zqs(ngrid)
    69         real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs
     63        DOUBLE PRECISION zqs(ngrid), zqs_tmp(ngrid), zt_temp(ngrid)
     64        real Lcp(ngrid,nlayer)
     65        real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,zqs_temp,zdqs
    7066        real zqs_temp_1, zqs_temp_2
    7167        integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice of generic_tracer
    … …  
    7773        DOUBLE PRECISION zx_q(ngrid)
    7874        DOUBLE PRECISION zy_q(ngrid)
    79         LOGICAL,SAVE :: firstcall=.true.
    80 !$OMP THREADPRIVATE(firstcall)
    81         IF (firstcall) THEN
    82                 write(*,*) "value for metallicity? "
    83                 metallicity=0.0 ! default value
    84                 call getin_p("metallicity",metallicity)
    85                 write(*,*) " metallicity = ",metallicity
    86 
    87                 write(*,*) "Deep generic tracer vapor mixing ratio ? (no effect if negative) "
    88                 qvap_deep=-1. ! default value
    89                 call getin_p("qvap_deep",qvap_deep)
    90                 write(*,*) " qvap_deep = ",qvap_deep   
    91 
    92                 write(*,*) "top generic tracer vapor mixing ratio ? (no effect if negative) "
    93                 qvap_top=-1. ! default value
    94                 call getin_p("qvap_top",qvap_top)
    95                 write(*,*) " qvap_top = ",qvap_top
    96 
    97                 write(*,*) " align_strato_cold_trap ? "
    98                 align_strato_cold_trap=.false. ! default value
    99                 call getin_p("align_strato_cold_trap",align_strato_cold_trap)
    100                 write(*,*) " align_strato_cold_trap = ",align_strato_cold_trap
    101 
    102                 firstcall = .false.
    103         ENDIF
     75        DOUBLE PRECISION zzx_q(ngrid)
     76        DOUBLE PRECISION zzy_q(ngrid)
     77        DOUBLE PRECISION zzt(ngrid)
     78
    10479!       Initialisation of outputs and local variables
    10580        pdtlsc(1:ngrid,1:nlayer)  = 0.0
    … …  
    125100                        metallicity_coeff=constants_metallicity_coeff(iq)
    126101
    127                         Lcp=RLVTT_generic/cpp ! need to be init here
     102                        Lcp(:,:)=RLVTT_generic/cpp(:,:) ! need to be init here
    128103
    129104                        !  Vertical loop (from top to bottom)
    130105                        DO k = nlayer, 1, -1
    131106                                zt(1:ngrid)=pt(1:ngrid,k)+pdt(1:ngrid,k)*ptimestep
     107                                zzt(1:ngrid)=pt(1:ngrid,k)+pdt(1:ngrid,k)*ptimestep
    132108
    133109                                ! Computes Psat and the partial condensation
    … …  
    139115                                        endif
    140116                                        zx_q(i) = pq(i,k,igcm_generic_vap)+pdq(i,k,igcm_generic_vap)*ptimestep
     117                                        zzx_q(i) = zx_q(i)
    141118
    142119                                        call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
    143120                                        zy_q(i) = pq(i,k,igcm_generic_ice)+pdq(i,k,igcm_generic_ice)*ptimestep
     121                                        zzy_q(i) = zy_q(i)
    144122                                       
    145123                                        if ((zx_q(i) .le. zqs_temp) .and. (zy_q(i) .eq. 0.)) then
    … …  
    153131
    154132                                                ! iterative process to stabilize the scheme when large water amounts JL12
     133                                               
     134                                                SELECT CASE (TRIM(thermo_phy))
     135                                               
     136                                                CASE ('thermo_uni_ideal')
     137                                               
    155138                                                zcond(i) = 0.0d0
    156139                                                Do nn=1,nitermax 
    157140                                                        call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
    158141                                                        zqs(i)=zqs_temp
    159                                                         call Lcpdqsat_generic(zt(i),local_p,psat_tmp,zqs_temp,zdqs,dlnpsat_tmp)
     142                                                        call Lcpdqsat_generic(zt(i),local_p,cpp(i,k),psat_tmp,zqs_temp,zdqs,dlnpsat_tmp)
    160143                                                        zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs)
    161144                                                        !zcond can be negative here
    162145                                                        zx_q(i) = zx_q(i) - zcond_iter
    163146                                                        zcond(i) = zcond(i) + zcond_iter
    164                                                         zt(i) = zt(i) + zcond_iter*Lcp
     147                                                        zt(i) = zt(i) + zcond_iter*Lcp(i,k)
    165148                                                        if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit
    166149                                                        if (nn.eq.nitermax) print*,'itermax in largescale'
    … …  
    173156                                                !                       so the line below is to check how much we can evaporate.
    174157                                                !                       If there is no ice available, zcond(i) will be 0. (first condidition of "if")
    175 
     158                                               
    176159                                                zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_generic_ice)+pdq(i,k,igcm_generic_ice)*ptimestep))
     160                                                zt(i) = zzt(i)+Lcp(i,k)*zcond(i)
     161                                                END SELECT
    177162                                       
    178163                                        endif
    … …  
    190175                                pdqvaplsc(1:ngrid,k,igcm_generic_vap)  = - zcond(1:ngrid)
    191176                                pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_vap)
    192                                 pdtlsc(1:ngrid,k)  = pdtlsc(1:ngrid,k) + pdqliqlsc(1:ngrid,k,igcm_generic_ice)*Lcp
     177                                pdtlsc(1:ngrid,k)  = pdtlsc(1:ngrid,k) + (zt(1:ngrid)-zzt(1:ngrid))/ptimestep
    193178
    194179                        Enddo ! k= nlayer, 1, -1

  • trunk/LMDZ.GENERIC/libf/phygeneric/condense_co2.F90

    r4077 r4146  
    1313      USE geometry_mod, only: latitude ! in radians
    1414      USE tracer_h, only: noms, rho_co2
    15       use comcstfi_mod, only: g, r, cpp
     15      use comcstfi_mod, only: g, rd_ref, cppd_ref
    1616
    1717      implicit none
    … …  
    159159         endif
    160160
    161          ccond=cpp/(g*latcond)
     161         ccond=cppd_ref/(g*latcond)
    162162         print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
    163163
    … …  
    287287               !w(ig,ilay,i_co2ice) = 0.0
    288288               w(ig,ilay,i_co2ice) = vstokes *  subptimestep * &
    289                    pplev(ig,ilay)/(r*pt(ig,ilay))
     289                   pplev(ig,ilay)/(rd_ref*pt(ig,ilay))
    290290
    291291            END DO

  • trunk/LMDZ.GENERIC/libf/phygeneric/conduction.F90

    r4033 r4146  
    88                        tsurf,zzlev,zzlay,muvar,qvar,firstcall,zdtconduc)
    99   
    10       use comcstfi_mod, only: r, cpp, mugaz
     10      use thermo_mod, only: cpp,r
    1111      use callkeys_mod, only: phitop_conduc,zztop,a_coeff,s_coeff,force_conduction
    1212      use conc_mod,     only: lambda
    … …  
    180180
    181181      DO i=1,nlayer-1
    182          muvol(:,i)=pplay(:,i)/(r*zt(:,i))
    183          alpha(:,i)=cpp*(muvol(:,i)/ptimestep) &
     182         muvol(:,i)=pplay(:,i)/(r(:,i)*zt(:,i))
     183         alpha(:,i)=cpp(:,i)*(muvol(:,i)/ptimestep) &
    184184                             *(zlev(:,i+1)-zlev(:,i))
    185185      ENDDO
    186186
    187       muvol(:,nlayer)=pplay(:,nlayer)/(r*zt(:,nlayer))
    188       alpha(:,nlayer)=cpp*(muvol(:,nlayer)/ptimestep) &
     187      muvol(:,nlayer)=pplay(:,nlayer)/(r(:,nlayer)*zt(:,nlayer))
     188      alpha(:,nlayer)=cpp(:,nlayer)*(muvol(:,nlayer)/ptimestep) &
    189189                             *(zlev(:,nlayer+1)-zlev(:,nlayer))
    190190

  • trunk/LMDZ.GENERIC/libf/phygeneric/convadj.F90

    r3893 r4146  
    77      subroutine convadj(ngrid,nlay,nq,ptimestep, &
    88                         pplay,pplev,ppopsk, &
    9                          pu,pv,ph,pq, &
    10                          pdufi,pdvfi,pdhfi,pdqfi, &
    11                          pduadj,pdvadj,pdhadj,pdqadj)
    12 
    13       use tracer_h, only: igcm_h2o_vap
    14       use comcstfi_mod, only: g
     9                         pu,pv,pt,ph,pq, &
     10                         pdufi,pdvfi,pdtfi,pdhfi,pdqfi, &
     11                         pduadj,pdvadj,pdtadj,pdhadj,pdqadj)
     12
     13      use tracer_h
     14      use comcstfi_mod, only: g, rd_ref, cppd_ref, cppv_ref
    1515      use generic_cloud_common_h, only: epsi_generic
    1616      use generic_tracer_index_mod, only: generic_tracer_index
    1717      use callkeys_mod, only: tracer,water,generic_condensation, &
    18                               virtual_theta_correction
     18                              virtual_theta_correction, thermo_phy
     19      use thermo_mod
    1920
    2021      implicit none
    … …  
    4243      REAL,INTENT(IN) :: ppopsk(ngrid,nlay) ! Exner (wrt surface pressure)
    4344      ! fields from dynamics
     45      REAL,INTENT(IN) :: pu(ngrid,nlay) ! zonal wind (m/s)
     46      REAL,INTENT(IN) :: pv(ngrid,nlay) ! meridional wind (m/s)
    4447      REAL,INTENT(IN) :: ph(ngrid,nlay) ! potential temperature (K)
    45       REAL,INTENT(IN) :: pu(ngrid,nlay) ! zonal wind (m/s)
    46       REAL,INTENT(IN) :: pv(ngrid,nlay) ! meridioanl wind (m/s)
     48      REAL,INTENT(IN) :: pt(ngrid,nlay) ! temperature (K)
    4749      REAL,INTENT(IN) :: pq(ngrid,nlay,nq) ! tracers (../kg_air)
    4850      ! tendencies due to previous physical processes
    … …  
    5052      REAL,INTENT(IN) :: pdvfi(ngrid,nlay) ! on meridional wind (m/s/s)
    5153      REAL,INTENT(IN) :: pdhfi(ngrid,nlay)! on potential temperature (/K/s)
     54      REAL,INTENT(IN) :: pdtfi(ngrid,nlay)! on temperature (/K/s)
    5255      REAL,INTENT(IN) :: pdqfi(ngrid,nlay,nq) ! on tracers (../kg_air/s)
    5356      ! tendencies due to convetive adjustement
    5457      REAL,INTENT(OUT) :: pduadj(ngrid,nlay) ! on zonal wind (m/s/s)
    5558      REAL,INTENT(OUT) :: pdvadj(ngrid,nlay) ! on meridinal wind (m/s/s)
     59      REAL,INTENT(OUT) :: pdtadj(ngrid,nlay) ! on temperature (/K/s)
    5660      REAL,INTENT(OUT) :: pdhadj(ngrid,nlay) ! on potential temperature (/K/s)
    5761      REAL,INTENT(OUT) :: pdqadj(ngrid,nlay,nq) ! on traceurs (../kg_air/s)
    … …  
    6165!     -----
    6266
    63       INTEGER ig,i,l,l1,l2,jj
    64       INTEGER jcnt, jadrs(ngrid)
    65 
    66       REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
    67       REAL zu(ngrid,nlay),zv(ngrid,nlay)
    68       REAL zh(ngrid,nlay), zvh(ngrid,nlay)
    69       REAL zu2(ngrid,nlay),zv2(ngrid,nlay)
    70       REAL zh2(ngrid,nlay),zhc(ngrid,nlay)
    71       REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
     67      INTEGER :: ig,i,l,l1,l2,jj,ll,k
     68      INTEGER :: jcnt, jadrs(ngrid)
     69
     70      DOUBLE PRECISION :: sig(ngrid,nlay+1),sdsig(ngrid,nlay),dsig(ngrid,nlay)
     71      REAL :: zu(ngrid,nlay),zv(ngrid,nlay)
     72      REAL :: zt(ngrid,nlay), ztv(ngrid,nlay)
     73      REAL :: zh(ngrid,nlay), zvh(ngrid,nlay)
     74      REAL :: zu2(ngrid,nlay),zv2(ngrid,nlay)
     75      REAL :: zt2(ngrid,nlay), ztv2(ngrid,nlay)
     76      REAL :: zh2(ngrid,nlay),zhc(ngrid,nlay)
     77      DOUBLE PRECISION :: zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
     78      REAL :: rloc(ngrid,nlay), rdloc(ngrid, nlay), rvloc(ngrid,nlay)
     79      REAL :: cpploc(ngrid,nlay), cppdloc(ngrid,nlay), cppvloc(ngrid,nlay)
     80      REAL :: rcploc(ngrid,nlay)
     81      REAL :: rvaj, rdaj, cppvaj, cppdaj, cppaj, rcpaj
    7282
    7383!     Tracers
    74       INTEGER iq
    75       REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
    76       REAL zqm(nq)
    77      
    78       integer igcm_generic_vap, igcm_generic_ice
    79       logical call_ice_vap_generic
    80 
    81       LOGICAL vtest(ngrid),down
     84      INTEGER :: iq
     85      REAL :: zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
     86      REAL :: zqm(nq), zqm2(nq)
     87     
     88      INTEGER :: igcm_generic_vap, igcm_generic_ice
     89      LOGICAL :: call_ice_vap_generic
     90
     91      LOGICAL :: vtest(ngrid),down
    8292
    8393!     for conservation test
    84       real masse,cadjncons
     94      REAL :: masse,cadjncons
    8595
    8696!     --------------
    … …  
    8898!     --------------
    8999      zh(:,:)=ph(:,:)+pdhfi(:,:)*ptimestep
     100      zt(:,:)=pt(:,:)+pdtfi(:,:)*ptimestep
    90101      zu(:,:)=pu(:,:)+pdufi(:,:)*ptimestep
    91102      zv(:,:)=pv(:,:)+pdvfi(:,:)*ptimestep
    … …  
    96107
    97108      zh2(:,:)=zh(:,:)
     109      zt2(:,:)=zt(:,:)
    98110      zu2(:,:)=zu(:,:)
    99111      zv2(:,:)=zv(:,:)
    100112      zq2(:,:,:)=zq(:,:,:)
     113
     114!----------------------------------------------------------------!
     115!---------beginning of select of the thermodynamics case---------!
     116!----------------------------------------------------------------!
     117     
     118      SELECT CASE (TRIM(thermo_phy))
     119     
     120      CASE ('thermo_uni_ideal')
    101121
    102122!     -----------------------------
    … …  
    151171!     Calculate sigma in this column
    152172          DO l=1,nlay+1
    153             sig(l)=pplev(i,l)/pplev(i,1)
     173            sig(i,l)=pplev(i,l)/pplev(i,1)
    154174       
    155175          ENDDO
    156176         DO l=1,nlay
    157             dsig(l)=sig(l)-sig(l+1)
    158             sdsig(l)=ppopsk(i,l)*dsig(l)
     177            dsig(i,l)=sig(i,l)-sig(i,l+1)
     178            sdsig(i,l)=ppopsk(i,l)*dsig(i,l)
    159179         ENDDO
    160180          l2 = 1
    … …  
    171191              l1 = l2 - 1
    172192              l  = l1
    173               zsm = sdsig(l2)
    174               zdsm = dsig(l2)
     193              zsm = sdsig(i,l2)
     194              zdsm = dsig(i,l2)
    175195              zhm = zh2(i, l2)
    176196              if (generic_condensation .and. virtual_theta_correction) then
    … …  
    181201
    182202              DO
    183                 zsm = zsm + sdsig(l)
    184                 zdsm = zdsm + dsig(l)
     203                zsm = zsm + sdsig(i,l)
     204                zdsm = zdsm + dsig(i,l)
    185205 
    186206                if (generic_condensation .and. virtual_theta_correction) then
    187                   zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
    188                   zhmc = zhmc + sdsig(l) * (zh2(i, l)*(1.e0+(1.e0/epsi_generic-1.e0)*zq(i,l,igcm_generic_vap)) - zhmc) / zsm
     207                  zhm = zhm + sdsig(i,l) * (zh2(i, l) - zhm) / zsm
     208                  zhmc = zhmc + sdsig(i,l) * (zh2(i, l)*(1.e0+(1.e0/epsi_generic-1.e0)*zq(i,l,igcm_generic_vap)) - zhmc) / zsm
    189209                else
    190                   zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
     210                  zhm = zhm + sdsig(i,l) * (zh2(i, l) - zhm) / zsm
    191211                  zhmc = zhm
    192212                endif
    … …  
    233253              end do
    234254              DO l = l1, l2
    235                 zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
     255                zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(i,l)
    236256                zh2(i, l) = zhm
    237257!     modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    238                 zum=zum+dsig(l)*zu2(i,l)
    239                 zvm=zvm+dsig(l)*zv2(i,l)
     258                zum=zum+dsig(i,l)*zu2(i,l)
     259                zvm=zvm+dsig(i,l)*zv2(i,l)
    240260!                zum=zum+dsig(l)*zu(i,l)
    241261!                zvm=zvm+dsig(l)*zv(i,l)
    242262                do iq=1,nq
    243                    zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq)
     263                   zqm(iq) = zqm(iq)+dsig(i,l)*zq2(i,l,iq)
    244264!                   zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq)
    245265!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    252272                end do
    253273              ENDDO
    254               zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
    255               zum=zum/(sig(l1)-sig(l2+1))
    256               zvm=zvm/(sig(l1)-sig(l2+1))
     274              zalpha=zalpha/(zhm*(sig(i,l1)-sig(i,l2+1)))
     275              zum=zum/(sig(i,l1)-sig(i,l2+1))
     276              zvm=zvm/(sig(i,l1)-sig(i,l2+1))
    257277              do iq=1,nq
    258                  zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1))
     278                 zqm(iq) = zqm(iq)/(sig(i,l1)-sig(i,l2+1))
    259279              end do
    260280
    … …  
    302322            print*,'cadjncons        = ',cadjncons
    303323         do l = 1, nlay
    304             print*,'dsig         = ',dsig(l)
     324            print*,'dsig         = ',dsig(i,l)
    305325         end do         
    306326         do l = 1, nlay
    307             print*,'dsig         = ',dsig(l)
    308          end do
    309          do l = 1, nlay
    310             print*,'sig         = ',sig(l)
     327            print*,'sig         = ',sig(i,l)
    311328         end do
    312329         do l = 1, nlay
    … …  
    315332         do l = 1, nlay+1
    316333            print*,'pplev(ig,:)         = ',pplev(i,l)
    317          end do
    318          do l = 1, nlay
    319             print*,'ph(ig,:)         = ',ph(i,l)
    320          end do
    321          do l = 1, nlay
    322             print*,'ph(ig,:)         = ',ph(i,l)
    323334         end do
    324335         do l = 1, nlay
    … …  
    346357
    347358      ENDDO
     359     
     360      END SELECT
     361     
     362!----------------------------------------------------------------!
     363!------------end of select of the thermodynamics case------------!
     364!----------------------------------------------------------------!
    348365
    349366      pdhadj(:,:)=(zh2(:,:)-zh(:,:))/ptimestep
    350367      pduadj(:,:)=(zu2(:,:)-zu(:,:))/ptimestep
    351368      pdvadj(:,:)=(zv2(:,:)-zv(:,:))/ptimestep
     369      pdtadj(:,:)=(zt2(:,:)-zt(:,:))/ptimestep
    352370
    353371      if(tracer) then

  • trunk/LMDZ.GENERIC/libf/phygeneric/dyn1d/calcenergy_kcm.F90

    r1403 r4146  
    44use params_h, only : Rc
    55use watercommon_h, only : mH2O
    6 use comcstfi_mod, only: g, mugaz
     6use comcstfi_mod, only: g, mugaz_ref
    77implicit none
    88
    … …  
    3333  double precision cp_neutral
    3434
    35   m_n = mugaz/1000.
     35  m_n = mugaz_ref/1000.
    3636  m_v = mH2O/1000.
    3737

  • trunk/LMDZ.GENERIC/libf/phygeneric/dyn1d/initracer_1D.F90

    r3893 r4146  
    77    use tracer_h, only: noms, mmol
    88    use datafile_mod, only: datadir
    9     use comcstfi_mod, only: mugaz
     9    use comcstfi_mod, only: mugaz_ref
    1010
    1111    implicit none
    … …  
    210210        if (dont_overwrite) cycle
    211211        do ilay=1,nlayer
    212             pq(ilay,iq) = pqx(ilay,iq)*mmol(iq)/mugaz
     212            pq(ilay,iq) = pqx(ilay,iq)*mmol(iq)/mugaz_ref
    213213        end do
    214214    end do

  • trunk/LMDZ.GENERIC/libf/phygeneric/dyn1d/kcm1d.F90

    r4077 r4146  
    355355!!! - physical frequency: nevermind, in inifis this is a simple print
    356356
    357   cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
    358   CALL inifis(1,llm,nq,0,86400.0,1,1.0,latitude,longitude,cell_area,rad,g,r,cpp)
     357  cppd_ref=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
     358  CALL inifis(1,llm,nq,0,86400.0,1,1.0,latitude,longitude,cell_area,rad,g,rd_ref,cppd_ref)
    359359
    360360  if(.not.global1d)then

  • trunk/LMDZ.GENERIC/libf/phygeneric/dyn1d/kcmprof_fn.F90

    r3893 r4146  
    44use watercommon_h, only : mH2O
    55use gases_h
    6 use comcstfi_mod, only: mugaz, cpp, g
     6use comcstfi_mod, only: mugaz_ref, cppd_ref, g
    77use callkeys_mod, only: co2cond
    88implicit none
    … …  
    8181  !-------------------------------
    8282  ! assign input variables
    83   m_n     = dble(mugaz/1000.)
    84   cp_n    = cpp
     83  m_n     = dble(mugaz_ref/1000.)
     84  cp_n    = cppd_ref
    8585  ! modify/generalise later??
    8686

  • trunk/LMDZ.GENERIC/libf/phygeneric/dyn1d/rcm1d.F

    r4112 r4146  
    1919     &         obliquit,nres,z0,coefvis,coefir,
    2020     &         timeperi,e_elips,p_elips
    21       use comcstfi_mod, only: pi, cpp, rad, g, r,
    22      &                        mugaz, rcp, omeg
     21      use comcstfi_mod, only: pi, cppd_ref, rad, g, rd_ref,
     22     &                        mugaz_ref, rcp_ref, omeg
    2323      use time_phylmdz_mod, only: daysec, dtphys, nday,
    2424     &                            diagfi_output_rate
    … …  
    651651!!! - physical constants: nevermind, things are done allright below
    652652!!! - physical frequency: nevermind, in inifis this is a simple print
    653       cpp=-9999. ! dummy init for inifis, will be rewrite later on
    654       r=-9999.   ! dummy init for inifis, will be rewrite later on
     653      cppd_ref=-9999. ! dummy init for inifis, will be rewrite later on
     654      rd_ref=-9999.   ! dummy init for inifis, will be rewrite later on
    655655      CALL inifis(1,llm,nq,day0,daysec,nday,dtphys,
    656      .            latitude,longitude,cell_area,rad,g,r,cpp)
     656     .            latitude,longitude,cell_area,rad,g,rd_ref,cppd_ref)
    657657
    658658      nsoil=nsoilmx
    … …  
    664664!!! We check everything is OK.
    665665      PRINT *,"CHECK"
    666       PRINT *,"--> mugaz = ",mugaz
    667       PRINT *,"--> cpp = ",cpp
    668       r = 8.314463E+0 * 1000.E+0 / mugaz
    669       rcp = r / cpp
     666      PRINT *,"--> mugaz_ref = ",mugaz_ref
     667      PRINT *,"--> cppd_ref = ",cppd_ref
     668      rd_ref = 8.314463E+0 * 1000.E+0 / mugaz_ref
     669      rcp_ref = rd_ref / cppd_ref
    670670
    671671c output spectrum?
    … …  
    10061006c  profil de temperature au premier appel
    10071007c  --------------------------------------
    1008       pks=psurf**rcp
     1008      pks=psurf**rcp_ref
    10091009
    10101010      if (restart) then
    … …  
    11821182!                 s(ilayer)=(play(ilayer)/psurf)**rcp
    11831183!              else
    1184           s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
     1184          s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp_ref
    11851185!              endif
    11861186              !s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
    1187           h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
     1187          h(ilayer)=cppd_ref*temp(ilayer)/(pks*s(ilayer))
    11881188       ENDDO
    11891189

  • trunk/LMDZ.GENERIC/libf/phygeneric/evap_generic.F90

    r2725 r4146  
    33    Use generic_cloud_common_h
    44    Use tracer_h
     5    use comcstfi_mod, only: cppd_ref
    56    implicit none
    67!==================================================================
    … …  
    3839
    3940    !   Evaporate all the ice from the tracer
    40     zlvdcp = RLVTT_generic/cpp ! RLVTT_generic is the latent heat of vaporization (comes from generic_cloud_common_h attention)
     41    zlvdcp = RLVTT_generic/cppd_ref ! RLVTT_generic is the latent heat of vaporization (comes from generic_cloud_common_h attention)
    4142 
    4243    DO l=1,nlayer

  • trunk/LMDZ.GENERIC/libf/phygeneric/generic_cloud_common_h.F90

    r3768 r4146  
    1313    !============================================================================
    1414
    15     use comcstfi_mod, only: r, cpp, mugaz
     15    use comcstfi_mod, only: rd_ref, mugaz_ref
    1616    implicit none
    1717    real, save :: m  ! molecular mass of the specie (g/mol)
    … …  
    107107            print*,"Unknow species (not in Mg, Fe, Na, KCl, Cr, Mn or Zn)"
    108108        endif
    109         epsi_generic = m/mugaz
     109        epsi_generic = m/mugaz_ref
    110110    end subroutine specie_parameters
    111111
    … …  
    182182        write(*,*) 'RLVTT_generic', RLVTT_generic
    183183
    184         epsi_generic = m/mugaz
     184        epsi_generic = m/mugaz_ref
    185185
    186186    end subroutine specie_parameters_table
    … …  
    203203        real, intent(out) :: psat,qsat !saturation vapor pressure (Pa) and mass mixing ratio at saturation (kg/kg) of the layer
    204204
    205         psat = pref*exp(-delta_vapH/(r*mugaz/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa)
     205        psat = pref*exp(-delta_vapH/(rd_ref*mugaz_ref/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa)
    206206
    207207        if (psat .gt. p) then
    … …  
    212212    end subroutine Psat_generic
    213213
    214     subroutine Lcpdqsat_generic(T,p,psat,qsat,dqsat,dlnpsat)
     214    subroutine Lcpdqsat_generic(T,p,cpp,psat,qsat,dqsat,dlnpsat)
    215215        implicit none
    216216
    … …  
    226226        real, intent(in) :: T ! Temperature (K)
    227227        real, intent(in) :: p ! Pressure (Pa)
     228        real, intent(in) :: cpp ! Heat capacity (J/kg)
    228229        real, intent(in) :: psat ! Saturation vapor pressure (Pa)
    229230        real, intent(in) :: qsat ! Mass mixing ratio at saturation (kg/kg)
    … …  
    235236        real :: dummy ! used to store d(ln(psat))/dT
    236237
    237         dummy = delta_vapH/((r*mugaz/1000.)*(T**2))
     238        dummy = delta_vapH/((rd_ref*mugaz_ref/1000.)*(T**2))
    238239
    239240        if (psat .gt. p) then
    … …  
    241242        else
    242243            dqsat = (RLVTT_generic/cpp) *qsat*(p/(p-(1-epsi_generic)*psat))*dummy
     244            dlnpsat = (RLVTT_generic/cpp) * dummy
    243245        endif
    244 
    245         dlnpsat = (RLVTT_generic/cpp) * dummy
    246246
    247247    end subroutine Lcpdqsat_generic
    … …  
    271271
    272272
    273         Tsat = 1/(1/Tref - (r*mugaz/1000.)/delta_vapH*Log(p/Pref))
     273        Tsat = 1/(1/Tref - (rd_ref*mugaz_ref/1000.)/delta_vapH*Log(p/Pref))
    274274
    275275    end subroutine Tsat_generic

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

    r4130 r4146  
    1919  use time_phylmdz_mod, only: diagfi_output_rate, slow_diagfi, &
    2020                              init_time, daysec, dtphys
    21   use comcstfi_mod, only: rad, cpp, g, r, rcp, &
    22                           mugaz, pi, avocado
     21  use comcstfi_mod, only: rad, cppd_ref, cppv_ref, cppc_ref, g, rd_ref, rcp_ref, &
     22                          mugaz_ref, pi, avocado
    2323  use planete_mod, only: nres
    2424  use planetwide_mod, only: planetwide_sumval
    … …  
    9090  ! initialize constants in comcstfi_mod
    9191  rad=prad
    92   cpp=pcpp
     92  cppd_ref=pcpp
    9393  g=pg
    94   r=pr
    95   rcp=r/cpp
    96   mugaz=8.314463*1000./pr ! dummy init
     94  rd_ref=pr
     95  rcp_ref=rd_ref/cppd_ref
     96  mugaz_ref=8.314463*1000./pr ! dummy init
    9797  pi=2.*asin(1.)
    9898  avocado = 6.022141e23   ! added by RW
    … …  
    11031103     if (is_master) write(*,*)trim(rname)//": generic_condensation = ",generic_condensation
    11041104     
     1105     if (generic_condensation) then
     1106                if (is_master) write(*,*)trim(rname)//": value for metallicity? "
     1107                metallicity=0.0 ! default value
     1108                call getin_p("metallicity",metallicity)
     1109                if (is_master) write(*,*)trim(rname)//": metallicity = ",metallicity
     1110
     1111                if (is_master) write(*,*)trim(rname)//": Deep generic tracer vapor mixing ratio ? (no effect if negative) "
     1112                qvap_deep=-1. ! default value
     1113                call getin_p("qvap_deep",qvap_deep)
     1114                if (is_master) write(*,*)trim(rname)//": qvap_deep = ",qvap_deep   
     1115
     1116                if (is_master) write(*,*)trim(rname)//":top generic tracer vapor mixing ratio ? (no effect if negative) "
     1117                qvap_top=-1. ! default value
     1118                call getin_p("qvap_top",qvap_top)
     1119                if (is_master) write(*,*)trim(rname)//": qvap_top = ",qvap_top
     1120
     1121                if (is_master) write(*,*)trim(rname)//": align_strato_cold_trap ? "
     1122                align_strato_cold_trap=.false. ! default value
     1123                call getin_p("align_strato_cold_trap",align_strato_cold_trap)
     1124                if (is_master) write(*,*)trim(rname)//": align_strato_cold_trap = ",align_strato_cold_trap
     1125     endif
     1126     
    11051127     if (is_master) write(*,*)trim(rname)//": Generic rain of tracers ?"
    11061128     generic_rain=.false. !default value
    11071129     call getin_p("generic_rain",generic_rain)
    11081130     if (is_master) write(*,*)trim(rname)//": generic_rain = ",generic_rain
     1131     
     1132     if (generic_rain) then
     1133            if (is_master) write(*,*)trim(rname)//": re-evaporate precipitations?"
     1134            evap_prec_generic=.true. ! default value
     1135            call getin_p("evap_prec_generic",evap_prec_generic)
     1136            if (is_master) write(*,*)trim(rname)//": evap_prec_generic = ",evap_prec_generic
     1137
     1138            if (evap_prec_generic) then
     1139               if (is_master) write(*,*)trim(rname)//": multiplicative constant in reevaporation"
     1140               evap_coeff_generic=1.   ! default value
     1141               call getin_p("evap_coeff_generic",evap_coeff_generic)
     1142               if (is_master) write(*,*)trim(rname)//": evap_coeff_generic = ",evap_coeff_generic       
     1143            end if
     1144
     1145            if (is_master) write(*,*)trim(rname)//": Precipitation scheme to use?"
     1146            precip_scheme_generic=1 ! default value
     1147            call getin_p("precip_scheme_generic",precip_scheme_generic)
     1148            if (is_master) write(*,*)trim(rname)//": precip_scheme_generic = ",precip_scheme_generic
     1149
     1150            if (precip_scheme_generic.eq.1) then
     1151               if (is_master) write(*,*)trim(rname)//": rainthreshold_generic in simple scheme?"
     1152               rainthreshold_generic=0. ! default value
     1153               call getin_p("rainthreshold_generic",rainthreshold_generic)
     1154               if (is_master) write(*,*)trim(rname)//": rainthreshold_generic = ",rainthreshold_generic
     1155               
     1156            else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then
     1157               
     1158               if (is_master) write(*,*)trim(rname)//": cloud GCS saturation level in non simple scheme?"
     1159               cloud_sat_generic=2.6e-4   ! default value
     1160               call getin_p("cloud_sat_generic",cloud_sat_generic)
     1161               if (is_master) write(*,*)trim(rname)//": cloud_sat_generic = ",cloud_sat_generic
     1162           
     1163               if (is_master) write(*,*)trim(rname)//": precipitation timescale in non simple scheme?"
     1164               precip_timescale_generic=3600.  ! default value
     1165               call getin_p("precip_timescale_generic",precip_timescale_generic)
     1166               if (is_master) write(*,*)trim(rname)//": precip_timescale_generic = ",precip_timescale_generic
     1167
     1168            else if (precip_scheme_generic.eq.4) then
     1169               
     1170               if (is_master) write(*,*)trim(rname)//": multiplicative constant in Boucher 95 precip scheme"
     1171               Cboucher_generic=1.   ! default value
     1172               call getin_p("Cboucher_generic",Cboucher_generic)
     1173               if (is_master) write(*,*)trim(rname)//": Cboucher_generic = ",Cboucher_generic   
     1174
     1175            endif
     1176     endif
    11091177
    11101178     if (is_master) write(*,*)trim(rname)//": Include moist adjustement for GCS ?"
    … …  
    12071275       call getin_p("nvarlayer",nvarlayer)
    12081276     endif
    1209 
    1210      if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
    1211      if (is_master) write(*,*)'Predefined Mg from dynamics is ',mugaz,'amu'
    1212 
    1213      force_cpp=.false. ! default value
    1214      call getin_p("force_cpp",force_cpp)
    1215      if (force_cpp) then
    1216       if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp
    1217       if (is_master) write(*,*)trim(rname)//": force_cpp is deprecated.",&
    1218       "Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
    1219       call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
    1220      endif
    1221 
    1222      if (is_master) write(*,*)trim(rname)//&
    1223      ": where do you want your cpp/mugaz value to come from?",&
    1224      "=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
     1277     
     1278     if (is_master) write(*,*)trim(rname)//&
     1279     ": where do you want your cppd_ref/mugaz_ref value to come from?",&
     1280     "=> 0: dynamics (3d), 1: forced in callphys.def (1d), 2: computed from gases.def (1d)?"
    12251281     cpp_mugaz_mode = 0 ! default value
    12261282     call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
    … …  
    12351291        call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
    12361292     endif
    1237 
    1238      if (cpp_mugaz_mode == 1) then
    1239        mugaz = -99999.
     1293     
     1294     if (is_master) write(*,*)trim(rname)//&
     1295     ": CHOOSE THERMODYNAMICS MOD"
     1296     thermo_phy='thermo_uni_ideal' !default value
     1297     call getin_p("thermo_phy",thermo_phy)
     1298     if (is_master) write(*,*)trim(rname)//": thermo_phy = ",trim(thermo_phy)
     1299     if(thermo_phy.eq.'thermo_uni_ideal') then
     1300       if (cpp_mugaz_mode == 0) write(*,*) "cppd_ref and mugaz_ref fixed by the dynamics"
     1301       if (cpp_mugaz_mode == 1) then
     1302         write(*,*) "cppd_ref and mugaz_ref fixed by the physics (callphys.def)"
     1303         cppd_ref = -99999.
     1304         if (is_master) write(*,*)trim(rname)//&
     1305           ": DRY SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
     1306         call getin_p("cppd_ref",cppd_ref)
     1307         IF (cppd_ref.eq.-99999.) THEN
     1308             PRINT *, "cppd_ref must be set if thermo_phy = thermo_uni_ideal and cpp_mugaz_mode == 1"
     1309             STOP
     1310         ENDIF
     1311         mugaz_ref = -99999.
     1312         if (is_master) write(*,*)trim(rname)//&
     1313           ": MEAN MOLECULAR MASS in g mol-1 ?"
     1314         call getin_p("mugaz_ref",mugaz_ref)
     1315         IF (mugaz_ref.eq.-99999.) THEN
     1316           call abort_physic(rname,"mugaz_ref must be set if cpp_mugaz_mode = 1",1)
     1317         ENDIF
     1318       endif
     1319       if (cpp_mugaz_mode == 2) write(*,*) "cppd_ref and mugaz_ref calculated by calc_cpp_mugaz in physics"
     1320       call su_gases
     1321       call calc_cpp_mugaz
     1322       if (generic_condensation) then
     1323         cppv_ref = -99999.
     1324         if (is_master) write(*,*)trim(rname)//&
     1325           ": VAPOUR TRACER SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
     1326         call getin_p("cppv_ref",cppv_ref)
     1327         IF (cppv_ref.eq.-99999.) THEN
     1328             PRINT *, "cppv_ref must be set if generic_condensation = true"
     1329             STOP
     1330         ENDIF
     1331       endif
     1332     elseif(thermo_phy.eq.'thermo_binary_ideal') then
     1333       write(*,*) "WARNING: thermo_binary_ideal can be used only wit the DYNAMICO CORE or in 1D"
     1334       if (.not.tracer) then
     1335          call abort_physic(rname,"Error: We need a tracer for the thermo_binary_ideal mod!",1)
     1336       endif
     1337       if (cpp_mugaz_mode == 0) write(*,*) "cppd_ref and mugaz_ref fixed by the dynamics"
     1338       if (cpp_mugaz_mode == 1) then
     1339         write(*,*) "cppd_ref and mugaz_ref fixed by the physics (callphys.def)"
     1340         cppd_ref = -99999.
     1341         if (is_master) write(*,*)trim(rname)//&
     1342           ": DRY SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
     1343         call getin_p("cppd_ref",cppd_ref)
     1344         IF (cppd_ref.eq.-99999.) THEN
     1345             PRINT *, "cppd_ref must be set if thermo_phy = thermo_binary_ideal"
     1346             STOP
     1347         ENDIF
     1348         mugaz_ref = -99999.
     1349         if (is_master) write(*,*)trim(rname)//&
     1350           ": MEAN MOLECULAR MASS in g mol-1 ?"
     1351         call getin_p("mugaz_ref",mugaz_ref)
     1352         IF (mugaz_ref.eq.-99999.) THEN
     1353           call abort_physic(rname,"mugaz_ref must be set if cpp_mugaz_mode = 1",1)
     1354         ENDIF
     1355       endif
     1356       if (cpp_mugaz_mode == 2) write(*,*) "cppd_ref and mugaz_ref calculated by calc_cpp_mugaz in physics"
     1357       call su_gases
     1358       call calc_cpp_mugaz
     1359       cppv_ref = -99999.
    12401360       if (is_master) write(*,*)trim(rname)//&
    1241          ": MEAN MOLECULAR MASS in g mol-1 ?"
    1242        call getin_p("mugaz",mugaz)
    1243        IF (mugaz.eq.-99999.) THEN
    1244          call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
     1361         ": VAPOUR TRACER SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
     1362       call getin_p("cppv_ref",cppv_ref)
     1363       IF (cppv_ref.eq.-99999.) THEN
     1364           PRINT *, "cppv_ref must be set if thermo_phy = thermo_binary_ideal"
     1365           STOP
    12451366       ENDIF
    1246        cpp = -99999.
     1367       cppc_ref = -99999.
    12471368       if (is_master) write(*,*)trim(rname)//&
    1248          ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
    1249        call getin_p("cpp",cpp)
    1250        IF (cpp.eq.-99999.) THEN
    1251            call abort_physic(rname, "cpp must be set if cpp_mugaz_mode = 1", 1)
     1369         ": CONDENSED TRACER SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
     1370       call getin_p("cppc_ref",cppc_ref)
     1371       IF (cppc_ref.eq.-99999.) THEN
     1372           PRINT *, "cppc_ref must be set if thermo_phy = thermo_binary_ideal"
     1373           STOP
    12521374       ENDIF
    1253        if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
    1254        if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu'
    1255  
    1256      endif ! of if (cpp_mugaz_mode == 1)
    1257      call su_gases
    1258      call calc_cpp_mugaz
     1375     else
     1376       call abort_physic(rname,"thermodynamics mod not recognized",1)
     1377     endif
    12591378
    12601379     if (is_master) then

  • trunk/LMDZ.GENERIC/libf/phygeneric/iniwrite.F

    r3865 r4146  
    22
    33      use comsoil_h, only: mlayer, nsoilmx
    4       USE comcstfi_mod, only: g, mugaz, omeg, rad, rcp, pi
     4      USE comcstfi_mod, only: g, mugaz_ref, omeg, rad, rcp_ref, pi
    55      USE vertical_layers_mod, ONLY: ap,bp,aps,bps,pseudoalt
    66!      USE logic_mod, ONLY: fxyhypb,ysinus
    … …  
    376376      write(*,*)'iniwrite: nbp_lon,nbp_lat,nbp_lev,idayref',
    377377     & nbp_lon,nbp_lat,nbp_lev,idayref
    378       write(*,*)'iniwrite: rad,omeg,g,mugaz,rcp',
    379      & rad,omeg,g,mugaz,rcp
     378      write(*,*)'iniwrite: rad,omeg,g,mugaz_ref,rcp_ref',
     379     & rad,omeg,g,mugaz_ref,rcp_ref
    380380      write(*,*)'iniwrite: daysec,dtphys',daysec,dtphys
    381381

  • trunk/LMDZ.GENERIC/libf/phygeneric/iniwrite_specIR.F

    r1621 r4146  
    33      use radinc_h, only: L_NSPECTI
    44      use radcommon_h, only: WNOI,DWNI
    5       use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
     5      use comcstfi_mod, only: rad, omeg, g, mugaz_ref, rcp_ref, pi
    66      use time_phylmdz_mod, only: daysec, dtphys
    77!      USE logic_mod, ONLY: fxyhypb,ysinus
    … …  
    7373      tab_cntrl(6)  = omeg
    7474      tab_cntrl(7)  = g
    75       tab_cntrl(8)  = mugaz
    76       tab_cntrl(9)  = rcp
     75      tab_cntrl(8)  = mugaz_ref
     76      tab_cntrl(9)  = rcp_ref
    7777      tab_cntrl(10) = daysec
    7878      tab_cntrl(11) = dtphys

  • trunk/LMDZ.GENERIC/libf/phygeneric/iniwrite_specVI.F

    r1621 r4146  
    33      use radinc_h, only: L_NSPECTV
    44      use radcommon_h, only: WNOV,DWNV
    5       use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
     5      use comcstfi_mod, only: rad, omeg, g, mugaz_ref, rcp_ref, pi
    66      use time_phylmdz_mod, only: daysec, dtphys
    77!      USE logic_mod, ONLY: fxyhypb,ysinus
    … …  
    7373      tab_cntrl(6)  = omeg
    7474      tab_cntrl(7)  = g
    75       tab_cntrl(8)  = mugaz
    76       tab_cntrl(9)  = rcp
     75      tab_cntrl(8)  = mugaz_ref
     76      tab_cntrl(9)  = rcp_ref
    7777      tab_cntrl(10) = daysec
    7878      tab_cntrl(11) = dtphys

  • trunk/LMDZ.GENERIC/libf/phygeneric/largescale.F90

    r2871 r4146  
    77          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water,Lcpdqsat_water
    88      USE tracer_h
     9      use callkeys_mod, only: qvap_deep
    910      IMPLICIT none
    1011
    … …  
    4041!     Options du programme
    4142      REAL, SAVE :: ratqs   ! determine largeur de la distribution de vapeur
    42       REAL, SAVE :: qvap_deep   ! deep mixing ratio of water vapor when simulating bottom less planets
    43 !$OMP THREADPRIVATE(ratqs, qvap_deep)
     43!$OMP THREADPRIVATE(ratqs)
    4444
    4545!     Variables locales
    … …  
    7373         call getin_p("ratqs",ratqs)
    7474         write(*,*) " ratqs = ",ratqs
    75 
    76          write(*,*) "Deep water vapor mixing ratio ? (no effect if negative) "
    77          qvap_deep=-1. ! default value
    78          call getin_p("qvap_deep",qvap_deep)
    79          write(*,*) " qvap_deep = ",qvap_deep
    8075
    8176         firstcall = .false.

  • trunk/LMDZ.GENERIC/libf/phygeneric/mass_redistribution_mod.F90

    r3893 r4146  
    1313       USE planete_mod, only: bp
    1414       use comcstfi_mod, only: g
    15        USE callkeys_mod, ONLY: water
     15       USE callkeys_mod, ONLY: water, generic_condensation
     16       use generic_tracer_index_mod, only: generic_tracer_index
     17       use generic_cloud_common_h
    1618       
    1719       IMPLICIT NONE
    … …  
    9698      REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
    9799      REAL sigma(nlayer+1)
    98 
    99 !   local saved variables
    100       LOGICAL, SAVE :: firstcall=.true.
    101 !$OMP THREADPRIVATE(firstcall)
    102 
    103 !----------------------------------------------------------------------
    104 
    105 !   Initialisation
    106 !   --------------
    107 !
    108       IF (firstcall) THEN
    109          firstcall=.false.       
    110       ENDIF
    111 !
     100     
     101!    Generic condensation     
     102      integer igcm_generic_vap, igcm_generic_ice
     103      logical call_ice_vap_generic
     104
    112105!======================================================================
    113106!    Calcul of h2o condensation 
    … …  
    159152         enddo
    160153      endif
     154     
     155      if (generic_condensation) then
     156     
     157         do iq=1,nq
     158
     159           call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
     160
     161         end do ! do iq=1,nq loop on tracers
     162     
     163         do ig=1,ngrid
     164            call Tsat_generic(pplev(ig,1)+zdmass_sum(ig,1)*g*ptimestep,Tsat(ig))
     165         enddo
     166#ifndef MESOSCALE
     167         call writediagfi(ngrid,'Tsat','saturation temperature at surface','',2,Tsat)
     168#endif
     169         
     170         do ig=1,ngrid
     171            if (ztsrf(ig).gt.Tsat(ig)) then
     172               zmassboil(ig)=(ptsrf(ig)-Tsat(ig))*pcapcal(ig)/RLVTT_generic/ptimestep
     173               if ((zmassboil(ig)*ptimestep.gt.pqs(ig,igcm_generic_vap)).and.(nint(rnat(ig)).eq.1)) then
     174                  zmassboil(ig)=pqs(ig,igcm_generic_vap)/ptimestep
     175               endif
     176               zmassboil(ig)=zmassboil(ig)*0.0 !momentary, should be 1. JL12
     177               pdqsmr(ig,igcm_generic_vap)=-zmassboil(ig)
     178               pdtsrfmr(ig)=-zmassboil(ig)*RLVTT_generic/pcapcal(ig)
     179               ztsrf(ig)=ptsrf(ig)+pdtsrfmr(ig)*ptimestep
     180            else
     181               zmassboil(ig)=0.
     182               pdtsrfmr(ig)=0.
     183            endif
     184         enddo
     185      endif
    161186
    162187!     *************************
    … …  
    225250         zqm(1,1:nq)=0. ! most tracer do not condense !
    226251         if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface
     252         if (generic_condensation) zqm(1,igcm_generic_vap)=1. ! flux is 100% condensable at surface
    227253         
    228254!        Van Leer scheme:

  • trunk/LMDZ.GENERIC/libf/phygeneric/moistadj.F90

    r3769 r4146  
    33  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, RCPV, Psat_water, Lcpdqsat_water
    44  USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
    5   use comcstfi_mod, only: r
     5  use comcstfi_mod, only: rd_ref
    66
    77  implicit none
    … …  
    4444!      PARAMETER (t_coup=234.0)
    4545      REAL seuil_vap
    46       PARAMETER (seuil_vap=1.0E-14)
     46      PARAMETER (seuil_vap=1.0E-10)
    4747
    4848!     Local variables
    … …  
    136136            v_pratio  = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp)
    137137            v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
    138             gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
     138            gamcpdz(i,k) = ( (rd_ref/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
    139139                / (((1.- v_zqs) + v_zqs * RCPV/RCPD)*v_pratio + v_zqs  * v_dlnpsat)               
    140140         ENDDO
    … …  
    234234         v_pratio  = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp)
    235235         v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
    236          gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
     236         gamcpdz(i,k) = ( (rd_ref/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
    237237                / (((1.- v_zqs) + v_zqs * RCPV/RCPD)*v_pratio + v_zqs  * v_dlnpsat)               
    238238      ENDDO

  • trunk/LMDZ.GENERIC/libf/phygeneric/moistadj_generic.F90

    r3769 r4146  
    1 subroutine moistadj_generic(ngrid, nlayer, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
     1subroutine moistadj_generic(ngrid, nlayer, nq, pt, pdt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
    22
    33   use generic_cloud_common_h
    44   use generic_tracer_index_mod, only: generic_tracer_index
    55   use tracer_h
    6    use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity
    7    use comcstfi_mod, only: r, cpp, mugaz
    8    use callkeys_mod, only: moist_convection_inhibition
     6   use callkeys_mod, only: moist_convection_inhibition, thermo_phy, metallicity
     7   use thermo_mod
     8   use comcstfi_mod
    99
    1010   implicit none
    … …  
    2828
    2929      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! temperature (K)
     30      REAL,INTENT(IN) :: pdt(ngrid,nlayer)
    3031      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracer (kg/kg)
    3132      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
    … …  
    3738      REAL,INTENT(OUT) :: rneb(ngrid,nlayer) ! cloud fraction
    3839
    39       ! Options :
    40       !real, save :: metallicity ! metallicity of planet
    41       REAL,SAVE :: metallicity = 0.0
    42 !$OMP THREADPRIVATE(metallicity)
    43 
    4440      ! local variables
    4541      REAL zt(ngrid,nlayer)      ! temperature (K)
    4642      REAL zq(ngrid,nlayer)      ! humidite specifique (kg/kg)
     43      REAL zql(ngrid,nlayer)     ! condensates concentration (kg/kg)
     44      REAL zqnc(ngrid,nlayer,nq) ! other tracers concentration (kg/kg)
    4745
    4846      REAL d_t(ngrid,nlayer)     ! temperature increment
    4947      REAL d_q(ngrid,nlayer)     ! incrementation pour vapeur d'eau
    5048      REAL d_ql(ngrid,nlayer)    ! incrementation pour l'eau liquide
     49      REAL d_qnc(ngrid,nlayer,nq)! incrementation pour les autres traceurs
    5150
    5251      ! REAL t_coup
    … …  
    6564      LOGICAL itest(ngrid)
    6665      REAL delta_q(ngrid, nlayer)
     66      REAL cpploc(ngrid, nlayer)
    6767      DOUBLE PRECISION :: cp_new_t(nlayer), v_cptt(ngrid,nlayer)
    6868      REAL cp_delta_t(nlayer)
    6969      DOUBLE PRECISION :: v_cptj(nlayer), v_cptjk1, v_ssig
    70       REAL v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
     70      REAL v_p, v_t, v_zqs,v_zql,v_tt2,v_cptt2,v_pratio,v_dlnpsat
    7171      REAL zqs(ngrid,nlayer), zdqs(ngrid,nlayer),zpsat(ngrid,nlayer),zdlnpsat(ngrid,nlayer)
    7272      REAL zq1(ngrid), zq2(ngrid)
    73       DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayer)
     73      DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayer),gamadz(ngrid,2:nlayer)
    7474      DOUBLE PRECISION :: zdp, zdpm
     75      REAL :: pqnc(nq), local_qnc(ngrid,nlayer,nq), local_qlk1
    7576
    7677      real q_cri(ngrid,nlayer) ! moist convection inhibition criterion
    … …  
    7980      REAL zflo ! flotabilite
    8081
    81       DOUBLE PRECISION :: local_q(ngrid,nlayer),local_t(ngrid,nlayer)
    82 
    83       REAL zdelta, zcor, zcvm5
     82      DOUBLE PRECISION :: local_q(ngrid,nlayer),local_t(ngrid,nlayer),local_ql(ngrid,nlayer)
     83     
     84      DOUBLE PRECISION A, B, C, D, delta, x1, x2, temp, den
     85      double precision :: sumqs,sumqvqc,sumdqs,sumdqsa,sumdqst,sumden
     86      REAL :: sig(ngrid,nlayer+1),dsig(ngrid,nlayer)
    8487
    8588      REAL dEtot, dqtot, masse ! conservation diagnostics
    … …  
    8790
    8891!     Indices of generic vapour and ice tracers
    89       real,save :: RCPD=0.0
    9092      INTEGER,SAVE :: i_vap_generic=0  ! Generic Condensable Species vapour
    9193      INTEGER,SAVE :: i_ice_generic=0  ! Generic Condensable Species ice
    92 !$OMP THREADPRIVATE(i_vap_generic,i_ice_generic,RCPD)
     94!$OMP THREADPRIVATE(i_vap_generic,i_ice_generic)
    9395
    9496      LOGICAL,SAVE :: firstcall=.TRUE.
    … …  
    9698
    9799      IF (firstcall) THEN
    98 
    99          RCPD = cpp
    100 
    101          write(*,*) "value for metallicity? "
    102          metallicity=0.0 ! default value
    103          call getin_p("metallicity",metallicity)
    104          write(*,*) " metallicity = ",metallicity
    105100
    106101         do iq=1, nq
    … …  
    115110               Tref = constants_Tref(iq)
    116111
    117                write(*,*) noms(igcm_generic_vap),", q_cri at ", Tref, "K (in kg/kg): ", ( 1 / (1 - 1/epsi_generic)) * (r * mugaz/1000.) / delta_vapH * Tref
     112               write(*,*) noms(igcm_generic_vap),", q_cri at ", Tref, "K (in kg/kg): ", ( 1 / (1 - 1/epsi_generic)) * (rd_ref * mugaz_ref/1000.) / delta_vapH * Tref
    118113
    119114            endif
    … …  
    131126      ! GCM -----> subroutine variables
    132127      zq(1:ngrid,1:nlayer)    = pq(1:ngrid,1:nlayer,i_vap_generic)+ pdq(1:ngrid,1:nlayer,i_vap_generic)*ptimestep
    133       zt(1:ngrid,1:nlayer)    = pt(1:ngrid,1:nlayer)
     128      zql(1:ngrid,1:nlayer)   = pq(1:ngrid,1:nlayer,i_ice_generic)+ pdq(1:ngrid,1:nlayer,i_ice_generic)*ptimestep
     129      zqnc(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq)+ pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
     130      zqnc(1:ngrid,1:nlayer,i_vap_generic) = 0.0
     131      zqnc(1:ngrid,1:nlayer,i_ice_generic) = 0.0
     132      zt(1:ngrid,1:nlayer)    = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
    134133      pdqmana(1:ngrid,1:nlayer,1:nq)=0.0
    135134
    … …  
    139138               zq(i,k)=0.0
    140139            endif
    141          ENDDO
    142       ENDDO
    143      
    144       local_q(1:ngrid,1:nlayer) = zq(1:ngrid,1:nlayer)
    145       local_t(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
     140            if(zql(i,k).lt.0.)then
     141               zql(i,k)=0.0
     142            endif
     143         ENDDO
     144      ENDDO
     145     
     146      local_q(1:ngrid,1:nlayer)  = zq(1:ngrid,1:nlayer)
     147      local_ql(1:ngrid,1:nlayer) = zql(1:ngrid,1:nlayer)
     148      local_qnc(1:ngrid,1:nlayer,1:nq)= zqnc(1:ngrid,1:nlayer,1:nq)
     149      local_t(1:ngrid,1:nlayer)  = zt(1:ngrid,1:nlayer)
    146150      rneb(1:ngrid,1:nlayer) = 0.0
    147151      d_ql(1:ngrid,1:nlayer) = 0.0
    148152      d_t(1:ngrid,1:nlayer)  = 0.0
    149153      d_q(1:ngrid,1:nlayer)  = 0.0
     154     
     155      do k=1,nlayer
     156        sig(:,k) = pplev(:,k)/pplev(:,1)
     157      enddo
     158      sig(:,nlayer+1) = pplev(:,nlayer+1)/pplev(:,1)
     159     
     160      do k=1,nlayer
     161        dsig(:,k) = sig(:,k)-sig(:,k+1)
     162      enddo
    150163
    151164      !     Calculate v_cptt
    152165      DO k = 1, nlayer
    153166         DO i = 1, ngrid
    154             v_cptt(i,k) = RCPD * local_t(i,k)
     167            v_cptt(i,k) = cppd_ref * local_t(i,k)
    155168            v_t = MAX(local_t(i,k),15.)
    156169            v_p = pplay(i,k)
    157170
    158171            call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k))
    159             call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
    160 
    161          ENDDO
    162       ENDDO
    163 
    164       ! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
     172            call Lcpdqsat_generic(v_t,v_p,cppd_ref,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
     173
     174         ENDDO
     175      ENDDO
     176
     177!----------------------------------------------------------------!
     178!---------beginning of select of the thermodynamics case---------!
     179!----------------------------------------------------------------!
     180
     181      SELECT CASE (TRIM(thermo_phy))
     182     
     183      CASE ('thermo_uni_ideal')
     184     
     185!------------------------------------ Calculate Gamma * Cp * dz: (gamma is the critical gradient)
    165186      DO k = 2, nlayer
    166187         DO i = 1, ngrid
    167188            zdp = pplev(i,k)-pplev(i,k+1)
    168189            zdpm = pplev(i,k-1)-pplev(i,k)
    169             ! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
     190            ! gamcpdz(i,k) = ( ( rd_ref/cppd_ref /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
    170191            !     +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
    171192            !* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
    … …  
    176197            v_pratio  = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp)
    177198            v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
    178             gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT_generic*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
    179                / (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)*v_pratio + v_zqs  * v_dlnpsat)               
     199            gamcpdz(i,k) = ( (rd_ref/cppd_ref*v_cptt2*(1.- v_zqs) + RLVTT_generic*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
     200               / (((1.- v_zqs) + v_zqs * cppv_ref/cppd_ref)*v_pratio + v_zqs  * v_dlnpsat)               
    180201            ! Note that gamcpdz is defined as positive, so -gamcpdz is the real moist-adiabatic gradient [dT/dz]_ad
    181202         ENDDO
    … …  
    188209         DO k = 1, nlayer
    189210            DO i = 1, ngrid
    190                q_cri(i,k) = ( 1 / (1 - 1/epsi_generic)) * r * mugaz/1000. / delta_vapH * zt(i,k)
     211               q_cri(i,k) = ( 1 / (1 - 1/epsi_generic)) * (rd_ref*mugaz_ref/1000.) / delta_vapH * zt(i,k)
    191212            ENDDO
    192213         ENDDO
    … …  
    251272      ENDDO
    252273
     274      local_qlk1 = 0.0
     275      den=0.0
     276      pqnc(:) = 0.0
     277      do k = k1, k2
     278         local_qlk1 = local_qlk1 + ((local_q(i,k)+local_ql(i,k)) - zqs(i,k) - zdqs(i,k)*(v_cptjk1/v_ssig - v_cptj(k) - v_cptt(i,k))/RLVTT_generic)*dsig(i,k)
     279         den            = den + dsig(i,k) - (local_q(i,k)+local_ql(i,k))*dsig(i,k)
     280         do iq=1,nq
     281           pqnc(iq) = pqnc(iq) + zqnc(i,k,iq)*dsig(i,k)
     282         enddo
     283      enddo
     284      local_qlk1 = local_qlk1/den
    253285
    254286      ! this right here is where the adjustment is done???
    255287      DO k = k1, k2
     288         local_ql(i,k) = local_qlk1
     289         do iq=1,nq
     290           local_qnc(i,k,iq) = pqnc(iq)/(sig(i,k1)-sig(i,k2+1))
     291         enddo
    256292         cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
    257293         cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
    258294         v_cptt(i,k)=cp_new_t(k)
    259295         local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT_generic
    260          local_t(i,k) = cp_new_t(k) / RCPD
     296         local_t(i,k) = cp_new_t(k) / cppd_ref
    261297
    262298         v_t = MAX(local_t(i,k),15.)
    … …  
    264300
    265301         call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k))
    266          call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
     302         call Lcpdqsat_generic(v_t,v_p,cppd_ref,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
    267303
    268304      ENDDO
    … …  
    275311
    276312!      DO k = k1, k2
    277 !         v_cptt(i,k) = RCPD * local_t(i,k)
     313!         v_cptt(i,k) = cppd_ref * local_t(i,k)
    278314!         v_t = local_t(i,k)
    279315!         v_p = pplay(i,k)
    … …  
    286322         zdpm = pplev(i,k-1) - pplev(i,k)
    287323         zdp = pplev(i,k) - pplev(i,k+1)
    288          ! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
     324         ! gamcpdz(i,k) = ( ( rd_ref/cppd_ref /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
    289325         ! +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
    290326         ! * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
    … …  
    295331         v_pratio  = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp)
    296332         v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
    297          gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT_generic*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
    298                / (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)*v_pratio + v_zqs  * v_dlnpsat)               
     333         gamcpdz(i,k) = ( (rd_ref/cppd_ref*v_cptt2*(1.- v_zqs) + RLVTT_generic*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
     334               / (((1.- v_zqs) + v_zqs * cppv_ref/cppd_ref)*v_pratio + v_zqs  * v_dlnpsat)               
    299335      ENDDO
    300336
    … …  
    3333699999 CONTINUE ! loop over all the points
    334370
    335 !-----------------------------------------------------------------------
    336 ! Determine the cloud fraction (hypothesis: nebulosity occurs
    337 ! where GCS vapor is reduced by adjustment):
    338 
    339       DO k = 1, nlayer
    340       DO i = 1, ngrid
    341          IF (itest(i)) THEN
    342          delta_q(i,k) = local_q(i,k) - zq(i,k)
    343          IF (delta_q(i,k).LT.0.) rneb(i,k)  = 1.0
    344          ENDIF
    345       ENDDO
    346       ENDDO
    347 
    348 ! Distribute GCS condensates into cloudy liquid/solid condensates (hypothesis:
    349 ! liquid/solid condensates are distributed to areas where GCS vapor
    350 ! decreases and are distributed in proportion to this decrease):
    351       DO i = 1, ngrid
    352          IF (itest(i)) THEN
    353          zq1(i) = 0.0
    354          zq2(i) = 0.0
    355          ENDIF
    356       ENDDO
    357       DO k = 1, nlayer
    358       DO i = 1, ngrid
    359          IF (itest(i)) THEN
    360          zdp = pplev(i,k)-pplev(i,k+1)
    361          zq1(i) = zq1(i) - delta_q(i,k) * zdp
    362          zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
    363          ENDIF
    364       ENDDO
    365       ENDDO
    366       DO k = 1, nlayer
    367       DO i = 1, ngrid
    368          IF (itest(i)) THEN
    369             IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
    370          ENDIF
    371       ENDDO
    372       ENDDO
     371
     372     END SELECT
     373
     374!----------------------------------------------------------------!
     375!------------end of select of the thermodynamics case------------!
     376!----------------------------------------------------------------!
     377
    373378
    374379      DO k = 1, nlayer
    … …  
    378383      ENDDO
    379384
    380       DO k = 1, nlayer
    381       DO i = 1, ngrid
    382          d_t(i,k) = local_t(i,k) - zt(i,k)
    383          d_q(i,k) = local_q(i,k) - zq(i,k)
    384       ENDDO
    385       ENDDO
    386 
    387385      ! now subroutine -----> GCM variables
    388386      DO k = 1, nlayer
    389387         DO i = 1, ngrid
    390            
    391             pdtmana(i,k)       = d_t(i,k)/ptimestep
    392             pdqmana(i,k,i_vap_generic) = d_q(i,k)/ptimestep
    393             pdqmana(i,k,i_ice_generic) = d_ql(i,k)/ptimestep
     388            pdtmana(i,k)       = (local_t(i,k) - zt(i,k))/ptimestep
     389            pdqmana(i,k,:)     = (local_qnc(i,k,:) - zqnc(i,k,:))/ptimestep
     390            pdqmana(i,k,i_vap_generic) = (local_q(i,k) - zq(i,k))/ptimestep
     391            pdqmana(i,k,i_ice_generic) = (local_ql(i,k) - zql(i,k))/ptimestep
    394392         
    395393         ENDDO

  • trunk/LMDZ.GENERIC/libf/phygeneric/molvis.F90

    r4033 r4146  
    99                        pvel,zzlev,zzlay,zdvelmolvis)
    1010     
    11       use comcstfi_mod, only: r, cpp, mugaz
    1211      use callkeys_mod, only: phitop_molvis,zztop
    1312      use conc_mod,     only: lambda
    1413      use gases_h
     14      use thermo_mod
    1515
    1616!=======================================================================
    … …  
    5151
    5252      INTEGER l,ig, nz     ! layer index, grid index, and number of layers
    53       real fac             ! conversion factor between thermal conductivity and molecular viscosity
     53      real fac(ngrid,nlayer)  ! conversion factor between thermal conductivity and molecular viscosity
    5454      REAL zvel(nlayer)    ! wind
    5555      real zt(nlayer)      ! temperature (K)
    … …  
    7171      nz = nlayer
    7272
    73       fac = (9 * cpp - 5 * (cpp - r)) / 4
     73      fac(:,:) = (9 * cpp(:,:) - 5 * (cpp(:,:) - r(:,:))) / 4
    7474
    7575      do ig=1,ngrid
    … …  
    8989        zlev(nz+1) = zlev(nz) + zztop
    9090       
    91         mu(1) = lambda(ig,1) / fac
     91        mu(1) = lambda(ig,1) / fac(ig,1)
    9292
    9393        DO l=2,nz
    94           mu(l)=lambda(ig,l)/fac
     94          mu(l)=lambda(ig,l)/fac(ig,l)
    9595        ENDDO
    9696   
    9797        DO l=1,nz-1
    98           muvol(l) = pplay(ig,l) / (r*zt(l))
     98          muvol(l) = pplay(ig,l) / (r(ig,l)*zt(l))
    9999          alpha(l) = (muvol(l) / ptimestep) * (zlev(l+1) - zlev(l))
    100100        ENDDO
    101101
    102         muvol(nz) = pplay(ig,nz) / (r*zt(nz))
     102        muvol(nz) = pplay(ig,nz) / (r(ig,nz)*zt(nz))
    103103        alpha(nz) = (muvol(nz) / ptimestep) * (zlev(nz+1) - zlev(nz))
    104104

  • trunk/LMDZ.GENERIC/libf/phygeneric/newsedim.F

    r4079 r4146  
    99     
    1010      use ioipsl_getin_p_mod, only: getin_p
    11       use comcstfi_mod, only: r, g
     11      use comcstfi_mod, only: g, rd_ref
    1212      use gases_h, only: gfrac, gnom, ngasmx
    1313      use gases_h, only: igas_CH4, igas_CO2, igas_H2, igas_H2O
    … …  
    196196
    197197c Correction for high Reynolds number
    198             Reynolds=2. * pplev(ig,l) / r / pt(ig,l) *
     198            Reynolds=2. * pplev(ig,l) / rd_ref / pt(ig,l) *
    199199     &        rsurf * vstokes(ig,l) / visc(ig,l)
    200200            if(Reynolds.ge.1.0) then
    … …  
    202202              Cd=24. / Reynolds * (1. + 0.15 * Reynolds**0.687) +
    203203     &           0.42 / (1. + 42500 / Reynolds**1.16)                       ! (Formula from Bagheri 2018)
    204               vstokes(ig,l) =(8./3.*pplev(ig,l)/r/pt(ig,l)*g*rfall**3 /
     204              vstokes(ig,l) =(8./3.*pplev(ig,l)/rd_ref/pt(ig,l)*g*rfall**3 /
    205205     &           rsurf**2/rho/Cd *
    206206     &           (1.+1.333*(a*pt(ig,l)/pplev(ig,l))/rsurf))**0.5
    207               Reynolds=2. * pplev(ig,l) / r / pt(ig,l) *
     207              Reynolds=2. * pplev(ig,l) / rd_ref / pt(ig,l) *
    208208     &           rsurf * vstokes(ig,l) / visc(ig,l)
    209209              enddo
    … …  
    234234cc             write(*,*) 'OK simple method l,w =', l, w(ig,l)
    235235cc            write(*,*) 'OK simple method dztop =', dztop
    236            w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l)))
     236           w(ig,l) = 1. - exp(-dztop*g/(rd_ref*pt(ig,l)))
    237237             !!! Diagnostic: JF. Fix: AS. Date: 05/11
    238238             !!! Probleme arrondi avec la quantite ci-dessus
    … …  
    242242
    243243             IF ( w(ig,l) .eq. 0. ) THEN
    244                 w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g
     244                w(ig,l) = ( dztop*g/(rd_ref*pt(ig,l)) ) * pplev(ig,l) /g
    245245             ELSE
    246246                w(ig,l) = w(ig,l) * pplev(ig,l) / g
    … …  
    264264               Ep = Ep - epaisseur(ig,l+k)
    265265!             ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
    266              ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
     266             ptop=exp(-(dztop-Ep)*g/(rd_ref*pt(ig,l+k)))
    267267             IF ( ptop .eq. 1. ) THEN
    268268                !PRINT*, 'newsedim: exposant trop petit ', ig, l
    269                 ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
     269                ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/
     270     &                    (rd_ref*pt(ig,l+k)))
    270271             ELSE
    271272                ptop=pplev(ig,l+k) * ptop

  • trunk/LMDZ.GENERIC/libf/phygeneric/nonoro_gwd_ran_mod.F90

    r2595 r4146  
    4747
    4848
    49      use comcstfi_mod, only: g, pi, r
     49     use comcstfi_mod, only: g, pi, rd_ref
    5050     USE ioipsl_getin_p_mod, ONLY : getin_p
    5151     use vertical_layers_mod, only : presnivs
    … …  
    195195
    196196        ! Characteristic vertical scale height
    197         H0 = r * tr / g
     197        H0 = rd_ref * tr / g
    198198        ! Control
    199199        if (deltat .LT. dtime) THEN
    … …  
    249249     DO LL = 2, nlayer + 1
    250250       !ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
    251         H0bis(:, LL-1) = r * TT(:, LL-1) / g
     251        H0bis(:, LL-1) = rd_ref * TT(:, LL-1) / g
    252252          ZH(:, LL) = ZH(:, LL-1) &
    253253           + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)

  • trunk/LMDZ.GENERIC/libf/phygeneric/phyredem.F90

    r3552 r4146  
    2121  use planete_mod, only: year_day, periastr, apoastr, peri_day, &
    2222                         obliquit, z0
    23   use comcstfi_mod, only: rad, omeg, g, mugaz, rcp
     23  use comcstfi_mod, only: rad, omeg, g, mugaz_ref, rcp_ref
    2424  use time_phylmdz_mod, only: daysec
    2525
    … …  
    9898
    9999! Informations about the planet, used by dynamics and physics
    100   tab_cntrl(5) = rad      ! planet radius (m)
    101   tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
    102   tab_cntrl(7) = g        ! gravity (m.s-2)
    103   tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1)
    104   tab_cntrl(9) = rcp      !  = r/cp  (=kappa in the dynamics)
    105   tab_cntrl(10) = daysec  ! length of a solar day (s)
     100  tab_cntrl(5) = rad       ! planet radius (m)
     101  tab_cntrl(6) = omeg      ! rotation rate (rad.s-1)
     102  tab_cntrl(7) = g         ! gravity (m.s-2)
     103  tab_cntrl(8) = mugaz_ref ! Molar mass of the atmosphere (g.mol-1)
     104  tab_cntrl(9) = rcp_ref   !  = r/cp  (=kappa in the dynamics)
     105  tab_cntrl(10) = daysec   ! length of a solar day (s)
    106106
    107107  tab_cntrl(11) = phystep  ! physics time step (s)

  • trunk/LMDZ.GENERIC/libf/phygeneric/physiq_mod.F90

    r4127 r4146  
    6262      use planete_mod, only: apoastr, periastr, year_day, peri_day, &
    6363                            obliquit, nres, z0
    64       use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
     64      use comcstfi_mod, only: pi, g, rcp_ref, rd_ref, rad, mugaz_ref, cppd_ref
     65      use thermo_mod
    6566      use time_phylmdz_mod, only: daysec
    6667      use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
    … …  
    7879                              tracer, UseTurbDiff, water, watercond, &
    7980                              waterrain, generic_rain, global1d, szangle, &
    80                               moistadjustment, moistadjustment_generic, varspec, varspec_data, nvarlayer
     81                              moistadjustment, moistadjustment_generic, varspec, &
     82                              varspec_data, nvarlayer,thermo_phy, metallicity
    8183      use generic_tracer_index_mod, only: generic_tracer_index
    8284      use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
    … …  
    300302      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers.
    301303      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
     304      real rtlaymean                 ! Mean r x temperature for calculation of zzlev
    302305
    303306! VARIABLES for the thermal plume model
    … …  
    459462      real rneb_generic(ngrid,nlayer,nq) ! GCS cloud fraction (generic condensation).
    460463      real psat_tmp_generic
    461       real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
    462 !$OMP THREADPRIVATE(metallicity)
    463464
    464465      real reffrad_generic_zeros_for_wrf(ngrid,nlayer) !  !!! this is temporary, it is only a list of zeros, it will be replaced when a generic aerosol will be implemented
    … …  
    792793         endif
    793794
    794 !        Set metallicity for GCS
    795 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    796          metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
    797          call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
    798 
    799795!        Set some parameters for the thermal plume model
    800796!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    801797         if (calltherm) then
    802             CALL init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RW)
     798            CALL init_thermcell_mod(g, rcp_ref, rd_ref, pi, T_h2o_ice_liq, RW)
    803799         endif
    804800
    … …  
    849845      endif
    850846
     847      ! Initialization of the thermodynamics mod
     848      ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
     849      if ((generic_condensation)) then
     850        ! take into account generic condensable specie (GCS) effect on mean molecular weight
     851        do iq=1,nq
     852          call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
     853        enddo
     854      else
     855        igcm_generic_vap=1
     856      endif
     857      call thermodynamics(thermo_phy, pplay, pplev, pt, ngrid, nlayer, nq, pq, igcm_generic_vap, zh, zpopsk)
    851858
    852859! ------------------------------------------------------
    … …  
    914921      ! Compute geopotential between layers.
    915922      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    916       zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
    917       do l=1,nlayer
    918          zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
    919       enddo
    920 
    921       zzlev(1:ngrid,1)=0.
    922 
    923       do l=2,nlayer
    924          do ig=1,ngrid
     923     
     924      ! Calculation of zzlay & zzlev for the first layer 
     925
     926      DO ig=1,ngrid
     927         zzlay(ig,1)=-(log(pplay(ig,1)/pplev(ig,1)))*r(ig,1)*pt(ig,1)/glat(ig)
     928         zzlev(ig,1)=0.
     929         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
     930      ENDDO
     931
     932      ! Compute zzlay
     933      DO l=2,nlayer
     934        DO ig=1,ngrid
     935         ! compute "mean" R X "mean" temperature of the layer
     936           if((pt(ig,l) .eq. pt(ig,l-1)) .and. (r(ig,l) .eq. r(ig,l-1))) then
     937             rtlaymean=r(ig,l)*pt(ig,l)
     938           else
     939             rtlaymean=(r(ig,l)*pt(ig,l)- r(ig,l-1)*pt(ig,l-1))/log(r(ig,l)*pt(ig,l)/(r(ig,l-1)*pt(ig,l-1)))
     940           endif
     941           
     942           zzlay(ig,l)=zzlay(ig,l-1)-(log(pplay(ig,l)/pplay(ig,l-1))*rtlaymean/glat(ig))
     943        ENDDO
     944      ENDDO
     945     
     946      ! Compute zzlev
     947      DO l=2,nlayer
     948         DO ig=1,ngrid
    925949            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
    926950            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
    927951            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
    928          enddo
    929       enddo
     952         ENDDO
     953      ENDDO 
    930954
    931955      !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
    932 
     956     
    933957      zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
    934958
    935       ! Compute potential temperature
    936       ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
     959      ! Compute mass of cell area
    937960      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    938961      do l=1,nlayer
    939962         do ig=1,ngrid
    940             zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
    941             zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    942963            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
    943964            massarea(ig,l)=mass(ig,l)*cell_area(ig)
    … …  
    953974      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
    954975      do l=1,nlayer
    955          pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /  &
     976         pw(1:ngrid,l)=(pw(1:ngrid,l)*r(1:ngrid,l)*pt(1:ngrid,l)) /  &
    956977                       (pplay(1:ngrid,l)*cell_area(1:ngrid))
    957978      enddo
    … …  
    10271048               if(water) then
    10281049                  ! take into account water effect on mean molecular weight
    1029                   muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
    1030                   muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
     1050                  muvar(1:ngrid,1:nlayer)=mugaz_ref/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
     1051                  muvar(1:ngrid,nlayer+1)=mugaz_ref/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
    10311052               elseif(generic_condensation) then
    10321053                  ! take into account generic condensable specie (GCS) effect on mean molecular weight
    1033                   do iq=1,nq
    1034 
    1035                      call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    1036 
    1037                      if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
    1038 
    1039                         epsi_generic=constants_epsi_generic(iq)
    1040 
    1041                         muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
    1042                         muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
    1043 
    1044                      endif
    1045                   end do ! do iq=1,nq loop on tracers
     1054
     1055                    epsi_generic=constants_epsi_generic(igcm_generic_vap)
     1056
     1057                    muvar(1:ngrid,1:nlayer)=mugaz_ref/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
     1058                    muvar(1:ngrid,nlayer+1)=mugaz_ref/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
     1059
    10461060               elseif(varspec) then
    10471061                  !take into account fixed variable mean molecular weight
    … …  
    10551069
    10561070               else
    1057                   muvar(1:ngrid,1:nlayer+1)=mugaz
     1071                  muvar(1:ngrid,1:nlayer+1)=mugaz_ref
    10581072               endif
    10591073
    … …  
    12071221         fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
    12081222         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
     1223         
     1224         ! Update thermodynamics
     1225         call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    12091226
    12101227         ! Test of energy conservation
    12111228         !----------------------------
    12121229         if(enertest)then
    1213             call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
    1214             call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW)
     1230            call planetwide_sumval(cpp(:,:)*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
     1231            call planetwide_sumval(cpp(:,:)*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW)
    12151232            !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
    12161233            call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
    12171234            call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*cell_area(:)/totarea_planet,dEtotsLW)
    1218             dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
    1219             dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
     1235            dEzRadsw(:,:)=cpp(:,:)*mass(:,:)*zdtsw(:,:)
     1236            dEzRadlw(:,:)=cpp(:,:)*mass(:,:)*zdtlw(:,:)
    12201237            if (is_master) then
    12211238               print*,'---------------------------------------------------------------'
    … …  
    12461263                          ptimestep,capcal,                      &
    12471264                          pplay,pplev,zzlay,zzlev,z0,            &
    1248                           pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
     1265                          pu,pv,pt,zh,zpopsk,pq,tsurf,emis,qsurf,&
    12491266                          pdt,pdq,zflubid,                       &
    12501267                          zdudif,zdvdif,zdtdif,zdtsdif,          &
    … …  
    12591276                       ptimestep,capcal,lwrite,               &
    12601277                       pplay,pplev,zzlay,zzlev,z0,            &
    1261                        pu,pv,zh,pq,tsurf,emis,qsurf,          &
     1278                       pu,pv,pt,zh,pq,tsurf,emis,qsurf,          &
    12621279                       zdh,pdq,zflubid,                       &
    12631280                       zdudif,zdvdif,zdhdif,zdtsdif,          &
    … …  
    12871304         end if ! of if (tracer)
    12881305
    1289 
     1306         ! Update thermodynamics
     1307         call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
     1308         
    12901309         ! test energy conservation
    12911310         !-------------------------
    12921311         if(enertest)then
    12931312
    1294             dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
     1313            dEzdiff(:,:)=cpp(:,:)*mass(:,:)*zdtdif(:,:)
    12951314            do ig = 1, ngrid
    12961315               dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
    … …  
    13971416            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq)
    13981417         ENDIF
     1418         
     1419         ! Update thermodynamics
     1420         call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    13991421
    14001422      ENDIF ! end of 'calltherm'
    … …  
    14101432         zdvadj(1:ngrid,1:nlayer)=0.0
    14111433         zdhadj(1:ngrid,1:nlayer)=0.0
    1412 
     1434         zdtadj(1:ngrid,1:nlayer)=0.0
    14131435
    14141436         call convadj(ngrid,nlayer,nq,ptimestep,            &
    14151437                      pplay,pplev,zpopsk,                   &
    1416                       pu,pv,zh,pq,                          &
    1417                       pdu,pdv,zdh,pdq,                      &
    1418                       zduadj,zdvadj,zdhadj,                 &
    1419                       zdqadj)
     1438                      pu,pv,pt,zh,pq,                          &
     1439                      pdu,pdv,pdt,zdh,pdq,                  &
     1440                      zduadj,zdvadj,zdtadj,                 &
     1441                      zdhadj,zdqadj)
    14201442
    14211443         pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
    … …  
    14271449            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
    14281450         end if
     1451         
     1452         ! Update thermodynamics
     1453         call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    14291454
    14301455         ! Test energy conservation
    14311456         if(enertest)then
    1432             call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot)
     1457            call planetwide_sumval(cpp(:,:)*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot)
    14331458            if (is_master) print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
    14341459         endif
    … …  
    14781503         print*,'du_nonoro: max ', maxval(d_u_hin), 'min ', minval(d_u_hin)
    14791504         print*,'dv_nonoro: max ', maxval(d_v_hin), 'min ', minval(d_v_hin)
     1505         
     1506         ! Update thermodynamics
     1507         call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    14801508
    14811509      ENDIF ! of IF (calllott_nonoro)
    … …  
    15051533         dqsurf(1:ngrid,igcm_co2_ice) = dqsurf(1:ngrid,igcm_co2_ice) + zdqsurfc(1:ngrid)
    15061534
     1535         ! Update thermodynamics
     1536         call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    15071537
    15081538         ! test energy conservation
    15091539         if(enertest)then
    1510             call planetwide_sumval(cpp*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot)
     1540            call planetwide_sumval(cpp(:,:)*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot)
    15111541            call planetwide_sumval(capcal(:)*zdtsurfc(:)*cell_area(:)/totarea_planet,dEtots)
    15121542            if (is_master) then
    … …  
    15571587                                                     + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
    15581588                  pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
     1589                 
     1590                  ! Update thermodynamics
     1591                  call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    15591592
    15601593                  ! Test energy conservation.
    15611594                  if(enertest)then
    1562                      call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
     1595                     call planetwide_sumval(cpp(:,:)*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
    15631596                     call planetwide_maxval(dtmoist(:,:),dtmoist_max)
    15641597                     call planetwide_minval(dtmoist(:,:),dtmoist_min)
    1565                      madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:)
     1598                     madjdEz(:,:)=cpp(:,:)*mass(:,:)*dtmoist(:,:)
    15661599
    15671600                     do ig=1,ngrid
    1568                         madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
     1601                        madjdE(ig) = SUM(cpp(:,:)*mass(:,:)*dtmoist(:,:))
    15691602                     enddo
    15701603
    … …  
    15941627               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayer)
    15951628               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer)
     1629               
     1630               ! Update thermodynamics
     1631               call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    15961632
    15971633               ! Test energy conservation.
    15981634               if(enertest)then
    1599                   lscaledEz(:,:) = cpp*mass(:,:)*dtlscale(:,:)
     1635                  lscaledEz(:,:) = cpp(:,:)*mass(:,:)*dtlscale(:,:)
    16001636                  do ig=1,ngrid
    1601                      lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
     1637                     lscaledE(ig) = SUM(cpp(:,:)*mass(:,:)*dtlscale(:,:))
    16021638                  enddo
    1603                   call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot)
     1639                  call planetwide_sumval(cpp(:,:)*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot)
    16041640
    16051641                  if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
    … …  
    16411677               dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
    16421678               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
     1679               
     1680               ! Update thermodynamics
     1681               call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    16431682
    16441683               ! Test energy conservation.
    16451684               if(enertest)then
    16461685
    1647                   call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot)
     1686                  call planetwide_sumval(cpp(:,:)*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot)
    16481687                  if (is_master) print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
    1649                   call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)/totarea_planet*RLVTT/cpp,dItot_tmp)
    1650                   call planetwide_sumval(cell_area(:)*zdqssnow(:)/totarea_planet*RLVTT/cpp,dItot)
     1688                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)/totarea_planet*RLVTT/cpp(:,:),dItot_tmp)
     1689                  call planetwide_sumval(cell_area(:)*zdqssnow(:)/totarea_planet*RLVTT/cpp(:,1),dItot)
    16511690                  dItot = dItot + dItot_tmp
    16521691                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dVtot_tmp)
    1653                   call planetwide_sumval(cell_area(:)*zdqsrain(:)/totarea_planet*RLVTT/cpp,dVtot)
     1692                  call planetwide_sumval(cell_area(:)*zdqsrain(:)/totarea_planet*RLVTT/cpp(:,1),dVtot)
    16541693                  dVtot = dVtot + dVtot_tmp
    16551694                  dEtot = dItot + dVtot
    … …  
    17241763           ! increment values of temperature:
    17251764           pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtchim(1:ngrid,1:nlayer)
     1765           
     1766           ! Update thermodynamics
     1767           call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    17261768
    17271769         END IF  ! of IF (photochem)
    … …  
    17411783               ! Moist Convective Adjustment.
    17421784               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1743                call moistadj_generic(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
     1785               call moistadj_generic(ngrid,nlayer,nq,pt,pdt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
    17441786
    17451787               pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)     &
    17461788                                                  + dqmoist(1:ngrid,1:nlayer,1:nq)
    17471789               pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
     1790               
     1791               ! Update thermodynamics
     1792               call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    17481793
    17491794               ! Test energy conservation.
    17501795               if(enertest)then
    1751                   call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
     1796                  call planetwide_sumval(cpp(:,:)*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
    17521797                  call planetwide_maxval(dtmoist(:,:),dtmoist_max)
    17531798                  call planetwide_minval(dtmoist(:,:),dtmoist_min)
    1754                   madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:)
     1799                  madjdEz(:,:)=cpp(:,:)*mass(:,:)*dtmoist(:,:)
    17551800
    17561801                  do ig=1,ngrid
    1757                      madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
     1802                     madjdE(ig) = SUM(cpp(:,:)*mass(:,:)*dtmoist(:,:))
    17581803                  enddo
    17591804
    … …  
    17831828            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqvaplscale_generic(1:ngrid,1:nlayer,1:nq)
    17841829            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqcldlscale_generic(1:ngrid,1:nlayer,1:nq)
     1830           
     1831            ! Update thermodynamics
     1832            call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    17851833
    17861834            if(enertest)then
    17871835               do ig=1,ngrid
    1788                   genericconddE(ig) = cpp*SUM(mass(:,:)*dt_generic_condensation(:,:))
     1836                  genericconddE(ig) = SUM(cpp(:,:)*mass(:,:)*dt_generic_condensation(:,:))
    17891837               enddo
    17901838
    1791                call planetwide_sumval(cpp*massarea(:,:)*dt_generic_condensation(:,:)/totarea_planet,dEtot)
     1839               call planetwide_sumval(cpp(:,:)*massarea(:,:)*dt_generic_condensation(:,:)/totarea_planet,dEtot)
    17921840
    17931841               if (is_master) print*,'In generic condensation atmospheric energy change =',dEtot,' W m-2'
    … …  
    18421890
    18431891            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq)
     1892           
     1893            ! Update thermodynamics
     1894            call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    18441895
    18451896            ! Test energy conservation.
    18461897            if(enertest)then
    18471898
    1848                call planetwide_sumval(cpp*massarea(:,:)*zdtrain_generic(:,:)/totarea_planet,dEtot)
     1899               call planetwide_sumval(cpp(:,:)*massarea(:,:)*zdtrain_generic(:,:)/totarea_planet,dEtot)
    18491900               if (is_master) print*,'In rain_generic atmospheric T energy change       =',dEtot,' W m-2'
    18501901
    … …  
    18591910                  if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
    18601911
    1861                      call planetwide_sumval(massarea(:,:)*dq_rain_generic_cld(:,:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dItot_tmp)
    1862                      call planetwide_sumval(cell_area(:)*zdqssnow_generic(:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dItot)
     1912                     call planetwide_sumval(massarea(:,:)*dq_rain_generic_cld(:,:,igcm_generic_ice)/totarea_planet*RLVTT/cpp(:,:),dItot_tmp)
     1913                     call planetwide_sumval(cell_area(:)*zdqssnow_generic(:,igcm_generic_ice)/totarea_planet*RLVTT/cpp(:,1),dItot)
    18631914                     dItot = dItot + dItot_tmp
    1864                      call planetwide_sumval(massarea(:,:)*dq_rain_generic_vap(:,:,igcm_generic_vap)/totarea_planet*RLVTT/cpp,dVtot_tmp)
    1865                      call planetwide_sumval(cell_area(:)*zdqsrain_generic(:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dVtot)
     1915                     call planetwide_sumval(massarea(:,:)*dq_rain_generic_vap(:,:,igcm_generic_vap)/totarea_planet*RLVTT/cpp(:,:),dVtot_tmp)
     1916                     call planetwide_sumval(cell_area(:)*zdqsrain_generic(:,igcm_generic_ice)/totarea_planet*RLVTT/cpp(:,1),dVtot)
    18661917                     dVtot = dVtot + dVtot_tmp
    18671918                     dEtot = dItot + dVtot
    … …  
    19271978            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
    19281979            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
     1980           
     1981            ! Update thermodynamics
     1982            call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    19291983
    19301984            ! Test water conservation
    … …  
    19482002         ! Updating Atmospheric Mass and Tracers budgets.
    19492003         if(mass_redistrib) then
    1950 
    1951             zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) *     &
     2004         
     2005            if (water) then
     2006
     2007              zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) *     &
    19522008               (   zdqevap(1:ngrid,1:nlayer)                          &
    19532009                 + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
    19542010                 + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)             &
    19552011                 + dqvaplscale(1:ngrid,1:nlayer) )
    1956 
     2012                 
     2013            else
     2014           
     2015              zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) *     &
     2016               (   dq_rain_generic_vap(1:ngrid,1:nlayer,igcm_generic_vap)             &
     2017                 + dqmoist(1:ngrid,1:nlayer,igcm_generic_vap)             &
     2018                 + dqvaplscale_generic(1:ngrid,1:nlayer,igcm_generic_vap) )
     2019                 
     2020            endif
     2021           
    19572022            do ig = 1, ngrid
    19582023               zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
    … …  
    19752040            pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
    19762041            zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
     2042           
     2043            ! Update thermodynamics
     2044            call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    19772045
    19782046         endif
    … …  
    20922160         call radioactive_tracers(ngrid,nlayer,nq,ptimestep,pq,zdqradio)
    20932161         pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqradio(1:ngrid,1:nlayer,1:nq)
     2162         
     2163        ! Update thermodynamics
     2164        call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    20942165
    20952166      endif! end of if 'tracer'
    … …  
    21952266      ! Dynamical heating diagnostic.
    21962267      do ig=1,ngrid
    2197          fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp
     2268         fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:)*cpp(ig,:))
    21982269      enddo
    21992270
    … …  
    24132484          pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmolvis(1:ngrid,1:nlayer)
    24142485        endif
     2486       
     2487        ! Update thermodynamics
     2488        call thermodynamics(thermo_phy, pplay, pplev, pt+pdt*ptimestep, ngrid, nlayer, nq, pq+pdq*ptimestep, igcm_generic_vap,zh,zpopsk)
    24152489
    24162490      endif ! of if (callthermos)
    … …  
    25092583
    25102584            if (water) then
    2511                vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
     2585               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz_ref/mmol(igcm_h2o_vap)
    25122586               call wstats(ngrid,"vmr_h2ovapor",                             &
    25132587                           "H2O vapour volume mixing ratio","mol/mol",       &

  • trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk.F90

    r4079 r4146  
    3737      use tracer_h, only: igcm_h2o_ice, igcm_h2o_vap, igcm_co2_ice
    3838      use tracer_h, only: constants_epsi_generic
    39       use comcstfi_mod, only: pi, mugaz, cpp
     39      use comcstfi_mod, only: pi, mugaz_ref
     40      use thermo_mod, only: cpp
    4041      use callkeys_mod, only: varactive,diurnal,tracer,water,varfixed,satval, &
    4142                              diagdtau,kastprof,strictboundcorrk,specOLR, &
    4243                              CLFvarying,tplanckmin,tplanckmax,global1d, &
    43                               generic_condensation, aerovenus, nvarlayer, varspec
     44                              generic_condensation, aerovenus, nvarlayer, varspec, &
     45                              metallicity, qvap_deep
    4446      use rad_correlatedk_opacities_stellar_mod, &
    4547        only: rad_correlatedk_opacities_stellar
    … …  
    232234      integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
    233235      logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
    234       real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
    235 !$OMP THREADPRIVATE(metallicity)
    236       REAL, SAVE :: qvap_deep   ! deep mixing ratio of water vapor when simulating bottom less planets
    237 !$OMP THREADPRIVATE(qvap_deep)
    238236
    239237      REAL :: maxvalue,minvalue
    … …  
    884882         endif
    885883         
    886          ! Initial values equivalent to mugaz.
     884         ! Initial values equivalent to mugaz_ref.
    887885         DO l=1,nlayer
    888             muvarrad(2*l)   = mugaz
    889             muvarrad(2*l+1) = mugaz
     886            muvarrad(2*l)   = mugaz_ref
     887            muvarrad(2*l+1) = mugaz_ref
    890888         END DO
    891889
    … …  
    11381136! ----------------------------------------------------------------
    11391137
    1140          Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
     1138         Cmk= 0.01 * 1.0 / (glat(ig) * mugaz_ref * 1.672621e-27) ! q_main=1.0 assumed.
    11411139         glat_ig=glat(ig)
    11421140
    … …  
    12531251         DO l=2,L_NLAYRAD
    12541252            dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
    1255                 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
     1253                *glat(ig)/(cpp(ig,L_NLAYRAD+1-l)*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
    12561254            dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
    1257                 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
     1255                *glat(ig)/(cpp(ig,L_NLAYRAD+1-l)*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
    12581256         END DO     
    12591257
    12601258!     These are values at top of atmosphere
    12611259         dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
    1262              *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
     1260             *glat(ig)/(cpp(ig,L_NLAYRAD)*scalep*(plevrad(3)-plevrad(2)))
    12631261         dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
    1264              *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
     1262             *glat(ig)/(cpp(ig,L_NLAYRAD)*scalep*(plevrad(3)-plevrad(2)))
    12651263
    12661264      !  Optical thickness diagnostics (added by JVO)

  • trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_online_recombination.F90

    r4081 r4146  
    230230    SUBROUTINE rad_correlatedk_recombination_main(igrid,nlayer,pq,pplay,pt,qvar,tmid,pmid)
    231231 
    232       USE comcstfi_mod, only: mugaz
     232      USE comcstfi_mod, only: mugaz_ref
    233233      USE radinc_h
    234234      USE radcommon_h
    … …  
    324324         if (.not. found) then ! set pq to top value
    325325            do iq=1,nrecomb_qotf
    326               pqr(iq,ip) = pq(nlayer,otf2tot_idx(iq))*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
     326              pqr(iq,ip) = pq(nlayer,otf2tot_idx(iq))*mugaz_ref/mmol(otf2tot_idx(iq)) ! mol/mol
    327327            enddo
    328328         else
    329329            if (ilay==0) then ! set pq to bottom value
    330330               do iq=1,nrecomb_qotf
    331                  pqr(iq,ip) = pq(1,otf2tot_idx(iq))*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
     331                 pqr(iq,ip) = pq(1,otf2tot_idx(iq))*mugaz_ref/mmol(otf2tot_idx(iq)) ! mol/mol
    332332               enddo
    333333            else ! standard : interp pq between layers
    … …  
    335335               do iq=1,nrecomb_qotf
    336336                 pqr(iq,ip) = pq(ilay,otf2tot_idx(iq))**fact * pq(ilay+1,otf2tot_idx(iq))**(1.0-fact)
    337                  pqr(iq,ip) = pqr(iq,ip)*mugaz/mmol(otf2tot_idx(iq)) ! mol/mol
     337                 pqr(iq,ip) = pqr(iq,ip)*mugaz_ref/mmol(otf2tot_idx(iq)) ! mol/mol
    338338               enddo
    339339            endif ! if ilay==nlayer

  • trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_opacities_stellar.F90

    r4077 r4146  
    1313  use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2, &
    1414                     igas_CH4, igas_CO2, igas_O2
    15   use comcstfi_mod, only: g, r, mugaz
     15  use comcstfi_mod, only: g, rd_ref, mugaz_ref
    1616  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec, &
    1717                          rayleigh
    … …  
    135135     if(kastprof)then
    136136        dz(k) = dpr(k)*(1000.0d0*8.314463d0/muvar(k))*TMID(K)/(g*PMID(K))
    137         U(k)  = Cmk*DPR(k)*mugaz/muvar(k)
     137        U(k)  = Cmk*DPR(k)*mugaz_ref/muvar(k)
    138138     else
    139         dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
    140         U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in rad_correlatedk_opacities_thermal.F 
     139        dz(k) = dpr(k)*rd_ref*TMID(K)/(glat_ig*PMID(K))*mugaz_ref/muvar(k)
     140        U(k)  = Cmk*DPR(k)*mugaz_ref/muvar(k)     ! only Cmk line in rad_correlatedk_opacities_thermal.F 
    141141            !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
    142142     endif

  • trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_opacities_thermal.F90

    r4077 r4146  
    1414  use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, &
    1515                     igas_CH4, igas_CO2, igas_O2
    16   use comcstfi_mod, only: g, r, mugaz
     16  use comcstfi_mod, only: g, rd_ref, mugaz_ref
    1717  use callkeys_mod, only: kastprof,continuum,graybody,varspec
    1818  use rad_correlatedk_online_recombination_mod, only: corrk_recombin, gasi_recomb
    … …  
    147147     if(kastprof)then
    148148        dz(k) = dpr(k)*(1000.0d0*8.314463d0/muvar(k))*TMID(K)/(g*PMID(K))
    149         U(k)  = Cmk*DPR(k)*mugaz/muvar(k)
     149        U(k)  = Cmk*DPR(k)*mugaz_ref/muvar(k)
    150150     else
    151         dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
    152         U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in rad_correlatedk_opacities_thermal.F 
     151        dz(k) = dpr(k)*rd_ref*TMID(K)/(glat_ig*PMID(K))*mugaz_ref/muvar(k)
     152        U(k)  = Cmk*DPR(k)*mugaz_ref/muvar(k)     ! only Cmk line in rad_correlatedk_opacities_thermal.F 
    153153            !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
    154154     endif

  • trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_read_opacity_tables.F90

    r4083 r4146  
    3434      use radcommon_h, only : wrefvar,WNOI,WNOV
    3535      use datafile_mod, only: datadir
    36       use comcstfi_mod, only: mugaz
     36      use comcstfi_mod, only: mugaz_ref
    3737      use gases_h, only: gnom, ngasmx, &
    3838                         igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4, &
    … …  
    382382         call getin_p("kappa_VI",kappa_VI)
    383383         write(*,*)" kappa_VI = ",kappa_VI
    384          kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule         
     384         kappa_VI=kappa_VI*1.e4* mugaz_ref * 1.672621e-27        ! conversion from m^2/kg to cm^2/molecule         
    385385     
    386386!        constant absorption coefficient in IR
    … …  
    389389         call getin_p("kappa_IR",kappa_IR)
    390390         write(*,*)" kappa_IR = ",kappa_IR       
    391          kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule
     391         kappa_IR=kappa_IR*1.e4* mugaz_ref * 1.672621e-27        ! conversion from m^2/kg to cm^2/molecule
    392392
    393393         write(*,*)"graybody: Visible / Infrared separation set at band: IR=",nIR_limit,", VI=",nVI_limit

  • trunk/LMDZ.GENERIC/libf/phygeneric/rad_newton_cooling.F90

    r4077 r4146  
    11subroutine rad_netwon_cooling(ngrid,nlayer,mu0,sinlat,popsk,temp,pplay,pplev,dtrad,firstcall)
    22       
    3   use comcstfi_mod, only: rcp, pi
     3  use comcstfi_mod, only: rcp_ref, pi
    44  use callkeys_mod, only: tau_relax
    55  implicit none
    … …  
    8383        dT_EP  = 70.
    8484
    85         sig_trop=(T_trop/T_surf)**(1./rcp)
     85        sig_trop=(T_trop/T_surf)**(1./rcp_ref)
    8686
    8787        do l=1,nlayer

  • trunk/LMDZ.GENERIC/libf/phygeneric/rain.F90

    r4079 r4146  
    66  use aerosol_radius, only: aerosol_radius_h2o_liquid_ice_separate
    77  USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
    8   use comcstfi_mod, only: g, r
     8  use comcstfi_mod, only: g, rd_ref
    99  implicit none
    1010
    … …  
    215215                     zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass(i,k)/ptimestep !there was a bug here
    216216                     zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k))    & !default was 2.e-5
    217                         *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here
     217                        *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*rd_ref ! BC modif here
    218218                     zqevt = MAX (zqevt, 0.0)
    219219                     zqev  = MIN (zqev, zqevt)
    … …  
    268268               IF (rneb(i,k).GT.0.0) THEN
    269269                  zoliq(i) = ql(i,k)
    270                   zrho(i)  = pplay(i,k) / ( zt(i,k) * R )
     270                  zrho(i)  = pplay(i,k) / ( zt(i,k) * rd_ref )
    271271                  zdz(i)   = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
    272272                  zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)

  • trunk/LMDZ.GENERIC/libf/phygeneric/rain_generic.F90

    r4079 r4146  
    88   use aerosol_radius, only: aerosol_radius_h2o_liquid_ice_separate ! only used for precip_scheme_generic >=2
    99   use tracer_h
    10    use comcstfi_mod, only: g, r, cpp
     10   use comcstfi_mod, only: g, cppc_ref
     11   use thermo_mod
    1112   use generic_tracer_index_mod, only: generic_tracer_index
     13   use callkeys_mod, only: thermo_phy,metallicity,evap_prec_generic,evap_coeff_generic, &
     14                           precip_scheme_generic,rainthreshold_generic,cloud_sat_generic, &
     15                           precip_timescale_generic,Cboucher_generic
    1216   implicit none
    1317 
    … …  
    5155   real d_q(ngrid,nlayer)        ! GCS vapor increment
    5256   real d_ql(ngrid,nlayer)       ! liquid GCS / ice increment
     57   real qcloud(ngrid), qcloud_tmp(ngrid)
    5358
    5459   integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice fo GCS
    55 
    56    real, save :: RCPD ! equal to cpp
    57 
    58    real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
    59    !$OMP THREADPRIVATE(metallicity)
    6060
    6161   ! Subroutine options
    6262   real,parameter :: seuil_neb=0.001  ! Nebulosity threshold
    6363
    64    integer,save :: precip_scheme_generic      ! id number for precipitaion scheme
    65    
    66    ! for simple scheme  (precip_scheme_generic=1)
    67    real,save :: rainthreshold_generic                ! Precipitation threshold in simple scheme
    68    
    69    ! for sundquist scheme  (precip_scheme_generic=2-3)
    70    real,save :: cloud_sat_generic                    ! Precipitation threshold in non simple scheme
    71    real,save :: precip_timescale_generic             ! Precipitation timescale
    72    
    7364   ! for Boucher scheme  (precip_scheme_generic=4)
    74    real,save :: Cboucher_generic ! Precipitation constant in Boucher 95 scheme
    7565   real,parameter :: Kboucher=1.19E8
    7666   real,save :: c1
    7767   
    78    !$OMP THREADPRIVATE(precip_scheme_generic,rainthreshold_generic,cloud_sat_generic,precip_timescale_generic,Cboucher_generic,c1)
     68   !$OMP THREADPRIVATE(c1)
    7969
    8070   integer,parameter :: ninter=5 ! only used for precip_scheme_generic >=2
    81 
    82    logical,save :: evap_prec_generic ! Does the rain evaporate ?
    83    REAL,SAVE :: evap_coeff_generic ! multiplication evaporation constant. 1. gives the model of gregory et al.
    84    !$OMP THREADPRIVATE(evap_prec_generic,evap_coeff_generic)
    8571
    8672   ! for simple scheme : precip_scheme_generic=1
    … …  
    10288   real tnext(ngrid,nlayer)
    10389
    104    real dmass(ngrid,nlayer)   ! mass of air in each grid cell
     90   real dmassm2(ngrid,nlayer)   ! mass of air in each grid cell
    10591   real dWtot
    10692 
    … …  
    10995   integer, save :: i_vap_generic=0  ! GCS vapour
    11096   integer, save :: i_ice_generic=0  ! GCS ice
    111    !$OMP THREADPRIVATE(i_vap_generic,i_ice_generic)
    112 
    113    LOGICAL,save :: firstcall=.true.
    114    !$OMP THREADPRIVATE(firstcall)
     97!$OMP THREADPRIVATE(i_vap_generic,i_ice_generic)
    11598
    11699   ! to call only one time the ice/vap pair of a tracer
    … …  
    142125         RLVTT_generic=constants_RLVTT_generic(iq)
    143126
    144          RCPD = cpp
    145          
    146 
    147          if (firstcall) then
    148 
    149             metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
    150             call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
    151 
    152             i_vap_generic=igcm_generic_vap
    153             i_ice_generic=igcm_generic_ice
    154            
    155             write(*,*) "rain: i_ice_generic=",i_ice_generic
    156             write(*,*) "      i_vap_generic=",i_vap_generic
    157 
    158             write(*,*) "re-evaporate precipitations?"
    159             evap_prec_generic=.true. ! default value
    160             call getin_p("evap_prec_generic",evap_prec_generic)
    161             write(*,*) " evap_prec_generic = ",evap_prec_generic
    162 
    163             if (evap_prec_generic) then
    164                write(*,*) "multiplicative constant in reevaporation"
    165                evap_coeff_generic=1.   ! default value
    166                call getin_p("evap_coeff_generic",evap_coeff_generic)
    167                write(*,*) " evap_coeff_generic = ",evap_coeff_generic   
    168             end if
    169 
    170             write(*,*) "Precipitation scheme to use?"
    171             precip_scheme_generic=1 ! default value
    172             call getin_p("precip_scheme_generic",precip_scheme_generic)
    173             write(*,*) " precip_scheme_generic = ",precip_scheme_generic
    174 
    175             if (precip_scheme_generic.eq.1) then
    176                write(*,*) "rainthreshold_generic in simple scheme?"
    177                rainthreshold_generic=0. ! default value
    178                call getin_p("rainthreshold_generic",rainthreshold_generic)
    179                write(*,*) " rainthreshold_generic = ",rainthreshold_generic
    180            
    181             !else
    182             !   write(*,*) "precip_scheme_generic = ", precip_scheme_generic
    183             !   write(*,*) "this precip_scheme_generic has not been implemented yet,"
    184             !   write(*,*) "only precip_scheme_generic = 1 has been implemented"
    185                
    186             else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then
    187                
    188                write(*,*) "cloud GCS saturation level in non simple scheme?"
    189                cloud_sat_generic=2.6e-4   ! default value
    190                call getin_p("cloud_sat_generic",cloud_sat_generic)
    191                write(*,*) " cloud_sat_generic = ",cloud_sat_generic
    192            
    193                write(*,*) "precipitation timescale in non simple scheme?"
    194                precip_timescale_generic=3600.  ! default value
    195                call getin_p("precip_timescale_generic",precip_timescale_generic)
    196                write(*,*) " precip_timescale_generic = ",precip_timescale_generic
    197 
    198             else if (precip_scheme_generic.eq.4) then
    199                
    200                write(*,*) "multiplicative constant in Boucher 95 precip scheme"
    201                Cboucher_generic=1.   ! default value
    202                call getin_p("Cboucher_generic",Cboucher_generic)
    203                write(*,*) " Cboucher_generic = ",Cboucher_generic       
    204                
    205                c1=1.00*1.097/rhowater*Cboucher_generic*Kboucher 
    206 
    207             endif
    208 
    209             PRINT*, 'in rain_generic.F, ninter=', ninter
    210 
    211             firstcall = .false.
    212 
    213          endif ! of if (firstcall)
     127         i_vap_generic=igcm_generic_vap
     128         i_ice_generic=igcm_generic_ice
    214129
    215130         ! GCM -----> subroutine variables
    … …  
    220135               q(i,k)    = pq(i,k,i_vap_generic)+pdq(i,k,i_vap_generic)*ptimestep
    221136               ql(i,k)   = pq(i,k,i_ice_generic)+pdq(i,k,i_ice_generic)*ptimestep
    222 
    223137               !q(i,k)    = pq(i,k,i_vap_generic)!+pdq(i,k,i_vap_generic)
    224138               !ql(i,k)   = pq(i,k,i_ice_generic)!+pdq(i,k,i_ice_generic)
    … …  
    251165               call Psat_generic(zt(i,k),p_tmp,metallicity,psat_tmp,zqs(i,k)) ! calculates psat_tmp & zqs(i,k)
    252166               ! metallicity --- is not used here but necessary to call function Psat_generic
    253                call Lcpdqsat_generic(zt(i,k),p_tmp,psat_tmp,zqs(i,k),dqsat(i,k),dlnpsat(i,k)) ! calculates dqsat(i,k) & dlnpsat(i,k)
     167               call Lcpdqsat_generic(zt(i,k),p_tmp,cppd_ref,psat_tmp,zqs(i,k),dqsat(i,k),dlnpsat(i,k)) ! calculates dqsat(i,k) & dlnpsat(i,k)
    254168               call Tsat_generic(p_tmp,Tsat(i,k)) ! calculates Tsat(i,k)
    255169
    … …  
    260174         do k = 1, nlayer
    261175            do i = 1, ngrid
    262                dmass(i,k)=(pplev(i,k)-pplev(i,k+1))/g ! mass per m² in each layer
     176               dmassm2(i,k)=(pplev(i,k)-pplev(i,k+1))/g ! mass per m2 in each layer
    263177            enddo
    264178         enddo
    … …  
    275189
    276190                     if(zt(i,k).gt.Tsat(i,k))then ! if temperature of the layer box is greater than Tsat
    277                         ! treat the case where all liquid water should boil
    278                         zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*dmass(i,k)/RLVTT_generic/ptimestep,precip_rate(i)) ! on évapore
    279                         precip_rate(i)=MAX(precip_rate(i)-zqev,0.) ! we withdraw from precip_rate the evaporated part
    280                         d_q(i,k)=zqev/dmass(i,k)*ptimestep ! quantité évaporée
    281                         d_t(i,k) = - d_q(i,k) * RLVTT_generic/RCPD
     191                        ! treat the case where all condensables should boil
     192                        zqev=MIN((zt(i,k)-Tsat(i,k))*cppd_ref*dmassm2(i,k)/RLVTT_generic/ptimestep,precip_rate(i)) ! on évapore
     193                        precip_rate_tmp(i)= precip_rate(i) - zqev
     194                        precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0)
     195                        qcloud(i) = precip_rate(i)/dmassm2(i,k)*ptimestep
     196                        qcloud_tmp(i) = precip_rate_tmp(i)/dmassm2(i,k)*ptimestep
     197                        d_q(i,k)=zqev/dmassm2(i,k)*ptimestep ! quantité évaporée
     198                       
     199                        SELECT CASE (TRIM(thermo_phy))
     200                        CASE ('thermo_uni_ideal')
     201                        d_t(i,k) = d_t(i,k) - d_q(i,k) * RLVTT_generic/cppd_ref
     202                        END SELECT
     203                        precip_rate(i)  = precip_rate_tmp(i)
    282204                     else
     205                     
    283206                        ! zqev/zqevt are the maximum amount of vapour that we are allowed to add by evaporation of rain
    284                         zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass(i,k)/(ptimestep*(1.d0+dqsat(i,k)))
     207                        zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmassm2(i,k)/(ptimestep*(1.d0+dqsat(i,k)))
    285208                        !max evaporation to reach saturation implictly accounting for temperature reduction
    286209                        zqevt= MAX (0.0, evap_coeff_generic*2.0e-5*(1.0-q(i,k)/zqs(i,k))    & !default was 2.e-5
    287                               *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*R) ! BC modif here, is it R or r/(mu/1000) ?
     210                              *sqrt(precip_rate(i))*dmassm2(i,k)/pplay(i,k)*zt(i,k)*R(i,k)) ! BC modif here, is it R or r/(mu/1000) ?
    288211                        zqev  = MIN (zqev, zqevt)
    289212                        zqev  = MAX (zqev, 0.0) ! not necessary according to previous lines
    … …  
    292215                        precip_rate_tmp(i)= precip_rate(i) - zqev
    293216                        precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0)
    294 
    295                         d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmass(i,k)*ptimestep
    296                         d_t(i,k) = - d_q(i,k) * RLVTT_generic/RCPD
     217                        qcloud(i) = precip_rate(i)/dmassm2(i,k)*ptimestep
     218                        qcloud_tmp(i) = precip_rate_tmp(i)/dmassm2(i,k)*ptimestep
     219                       
     220                        d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmassm2(i,k)*ptimestep
     221                       
     222                        SELECT CASE (TRIM(thermo_phy))
     223                       
     224                        CASE ('thermo_uni_ideal')
     225                        d_t(i,k) = d_t(i,k) - d_q(i,k) * RLVTT_generic/cppd_ref
     226                        END SELECT
     227                       
    297228                        precip_rate(i)  = precip_rate_tmp(i)
     229                       
    298230                     end if
    299231#ifdef MESOSCALE
    300                      d_t(i,k) = d_t(i,k)+(pphi(i,k+1)-pphi(i,k))*precip_rate(i)*ptimestep/(RCPD*dmass(i,k))
     232                     d_t(i,k) = d_t(i,k)+(pphi(i,k+1)-pphi(i,k))*precip_rate(i)*ptimestep/(cppd_ref*dmassm2(i,k))
    301233                     ! JL22. Accounts for gravitational energy of falling precipitations (probably not to be used in the GCM
    302234                     !   where the counterpart is not included in the dynamics.)
    … …  
    321253                     if ((ql(i,k)/zneb(i)).gt.lconvert)then ! precipitate!
    322254                        d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
    323                         precip_rate(i)   = precip_rate(i) - d_ql(i,k)*dmass(i,k)/ptimestep
     255                        precip_rate(i)   = precip_rate(i) - d_ql(i,k)*dmassm2(i,k)/ptimestep
    324256                     endif
    325257                  endif
    … …  
    333265                  if (rneb(i,k).GT.0.0) then
    334266                     zoliq(i) = ql(i,k)
    335                      zrho(i)  = pplay(i,k) / ( zt(i,k) * R )
     267                     zrho(i)  = pplay(i,k) / ( zt(i,k) * R(i,k) )
    336268                     zdz(i)   = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
    337269                     zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
    … …  
    415347                  if (rneb(i,k).GT.0.0) then
    416348                     d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
    417                      precip_rate(i)   = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmass(i,k)/ptimestep
     349                     precip_rate(i)   = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmassm2(i,k)/ptimestep
    418350                  endif
    419351               enddo
    … …  
    449381               reevap_precip(i,i_vap_generic)=0.
    450382               do k=1,nlayer
    451                   reevap_precip(i,i_vap_generic)=reevap_precip(i,i_vap_generic)+dq_rain_generic_vap(i,k,i_vap_generic)*dmass(i,k)
     383                  reevap_precip(i,i_vap_generic)=reevap_precip(i,i_vap_generic)+dq_rain_generic_vap(i,k,i_vap_generic)*dmassm2(i,k)
    452384               enddo
    453385            enddo

  • trunk/LMDZ.GENERIC/libf/phygeneric/tabfi_mod.F90

    r3908 r4146  
    6363      use planete_mod, only: year_day, periastr, apoastr, peri_day, &
    6464                             obliquit, z0
    65       use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r
     65      use comcstfi_mod, only: rad, omeg, g, mugaz_ref, rcp_ref, cppd_ref, rd_ref
    6666      use time_phylmdz_mod, only: dtphys, daysec
    6767      use callkeys_mod, only: cpp_mugaz_mode
    … …  
    116116          call abort_physic(modname,"Missing value for omega in def files!",1)
    117117        endif
    118         mugaz=(8.314463/r)*1.E3
     118        mugaz_ref=(8.314463/rd_ref)*1.E3
    119119        ! daysec already set by inifis
    120120        ! dtphys alread set by inifis
    … …  
    189189      omeg = tab_cntrl(tab0+6)
    190190      g = tab_cntrl(tab0+7)
    191       mugaz = tab_cntrl(tab0+8)
    192       rcp = tab_cntrl(tab0+9)
    193       cpp=(8.314463/(mugaz/1000.0))/rcp
     191      mugaz_ref = tab_cntrl(tab0+8)
     192      rcp_ref = tab_cntrl(tab0+9)
     193      cppd_ref=(8.314463/(mugaz_ref/1000.0))/rcp_ref
    194194! daysec and dtphys are already set by inifis
    195195      cntrl_daysec = tab_cntrl(tab0+10)
    … …  
    235235      p_omeg = omeg
    236236      p_g = g
    237       p_cpp = cpp
    238       p_mugaz = mugaz
     237      p_cpp = cppd_ref
     238      p_mugaz = mugaz_ref
    239239      p_daysec = daysec
    240240      p_rad=rad
    … …  
    259259      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
    260260      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
    261       write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
    262       write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
     261      write(*,5) '(8)       mugaz_ref',tab_cntrl(tab0+8),mugaz_ref
     262      write(*,5) '(9)         rcp_ref',tab_cntrl(tab0+9),rcp_ref
    263263      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
    264264
    … …  
    316316      write(*,*) '(6)      omeg     : planet rotation rate (rad/s)'
    317317      write(*,*) '(7)      g        : gravity (m/s2)'
    318       write(*,*) '(8)      mugaz    : molecular mass '
     318      write(*,*) '(8)    mugaz_ref  : molecular mass '
    319319      write(*,*) '                       of the atmosphere (g/mol)'
    320       write(*,*) '(9)      rcp      : r/Cp'
     320      write(*,*) '(9)    rcp_ref    : r/Cp'
    321321      write(*,*) '(8)+(9)  calc_cpp_mugaz : r/Cp and mugaz '
    322322      write(*,*) '                 computed from gases.def'
    … …  
    477477          write(*,*) ' g (new value):',g
    478478
    479         else if (modif(1:len_trim(modif)).eq.'mugaz') then
    480           write(*,*) 'current value:',mugaz
    481           write(*,*) 'enter new value:'
    482  123      read(*,*,iostat=ierr) mugaz
     479        else if (modif(1:len_trim(modif)).eq.'mugaz_ref') then
     480          write(*,*) 'current value:',mugaz_ref
     481          write(*,*) 'enter new value:'
     482 123      read(*,*,iostat=ierr) mugaz_ref
    483483          if(ierr.ne.0) goto 123
    484484          write(*,*)
    485           write(*,*) ' mugaz (new value):',mugaz
    486           r=8.314463/(mugaz/1000.0)
    487           write(*,*) ' R (new value):',r
    488 
    489         else if (modif(1:len_trim(modif)).eq.'rcp') then
    490           write(*,*) 'current value:',rcp
    491           write(*,*) 'enter new value:'
    492  124      read(*,*,iostat=ierr) rcp
     485          write(*,*) ' mugaz_ref (new value):',mugaz_ref
     486          rd_ref=8.314463/(mugaz_ref/1000.0)
     487          write(*,*) ' rd_ref (new value):',rd_ref
     488
     489        else if (modif(1:len_trim(modif)).eq.'rcp_ref') then
     490          write(*,*) 'current value:',rcp_ref
     491          write(*,*) 'enter new value:'
     492 124      read(*,*,iostat=ierr) rcp_ref
    493493          if(ierr.ne.0) goto 124
    494494          write(*,*)
    495           write(*,*) ' rcp (new value):',rcp
    496           r=8.314463/(mugaz/1000.0)
    497           cpp=r/rcp
    498           write(*,*) ' cpp (new value):',cpp
     495          write(*,*) ' rcp_ref (new value):',rcp_ref
     496          rd_ref=8.314463/(mugaz_ref/1000.0)
     497          cppd_ref=rd_ref/rcp_ref
     498          write(*,*) ' cppd_ref (new value):',cppd_ref
    499499
    500500        else if (modif(1:len_trim(modif)).eq.'calc_cpp_mugaz') then
    501           write(*,*) 'current value rcp, mugaz:',rcp,mugaz
     501          write(*,*) 'current value rcp_ref, mugaz_ref:',rcp_ref,mugaz_ref
    502502          cpp_mugaz_mode = 2
    503503          call su_gases
    504504          call calc_cpp_mugaz
    505505          write(*,*)
    506           write(*,*) ' cpp (new value):',cpp
    507           write(*,*) ' mugaz (new value):',mugaz
    508           r=8.314463/(mugaz/1000.0)
    509           rcp=r/cpp
    510           write(*,*) ' rcp (new value):',rcp
     506          write(*,*) ' cppd_ref (new value):',cppd_ref
     507          write(*,*) ' mugaz_ref (new value):',mugaz_ref
     508          rd_ref=8.314463/(mugaz_ref/1000.0)
     509          rcp_ref=rd_ref/cppd_ref
     510          write(*,*) ' rcp_ref (new value):',rcp_ref
    511511         
    512512        else if (modif(1:len_trim(modif)).eq.'daysec') then
    … …  
    546546      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
    547547      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
    548       write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
    549       write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
     548      write(*,5) '(8)       mugaz_ref',tab_cntrl(tab0+8),mugaz_ref
     549      write(*,5) '(9)         rcp_ref',tab_cntrl(tab0+9),rcp_ref
    550550      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
    551551 
    … …  
    578578      p_omeg = omeg
    579579      p_g = g
    580       p_cpp = cpp
    581       p_mugaz = mugaz
     580      p_cpp = cppd_ref
     581      p_mugaz = mugaz_ref
    582582      p_daysec = daysec
    583583      p_rad=rad
    … …  
    589589      ! Initialize controle variable for XIOS & diagfi
    590590
    591         use comcstfi_mod,        only: g, mugaz, omeg, rad, rcp
     591        use comcstfi_mod,        only: g, mugaz_ref, omeg, rad, rcp_ref
    592592        use time_phylmdz_mod,  only: daysec, dtphys
    593593        use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
    … …  
    611611        tab_cntrl(6)  = omeg
    612612        tab_cntrl(7)  = g
    613         tab_cntrl(8)  = mugaz
    614         tab_cntrl(9)  = rcp
     613        tab_cntrl(8)  = mugaz_ref
     614        tab_cntrl(9)  = rcp_ref
    615615        tab_cntrl(10) = daysec
    616616        tab_cntrl(11) = dtphys

  • trunk/LMDZ.GENERIC/libf/phygeneric/thermcell_env.F90

    r3663 r4146  
    2626      USE watercommon_h, ONLY: RLvCp, RETV, Psat_water
    2727      USE tracer_h, ONLY: igcm_h2o_vap, igcm_h2o_ice
    28       USE callkeys_mod, ONLY: water, generic_condensation
    29       USE comcstfi_mod, ONLY: r, cpp, mugaz
     28      USE callkeys_mod, ONLY: water, generic_condensation, metallicity
     29      USE comcstfi_mod, ONLY: rd_ref, cppd_ref, mugaz_ref
    3030      USE generic_cloud_common_h, ONLY: Psat_generic, epsi_generic, RLVTT_generic
    3131      USE generic_tracer_index_mod, ONLY: generic_tracer_index
    … …  
    7575      INTEGER :: igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
    7676      LOGICAL :: call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
    77       REAL, SAVE :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
    78       !$OMP THREADPRIVATE(metallicity)
    7977      REAL :: RETV_generic, RV_generic, RLvCp_generic
    8078
    … …  
    9088      zqt(:,:) = 0.
    9189      zql(:,:) = 0.
    92      
    93       metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
    94        
     90         
    9591!===============================================================================
    9692! Condensation and latent heat release
    … …  
    132128               ENDDO
    133129
    134                RV_generic = (8.314463*1000.)/(epsi_generic*mugaz)
    135                RETV_generic = RV_generic/r-1.
    136                RLvCp_generic = RLVTT_generic/cpp
     130               RV_generic = (8.314463*1000.)/(epsi_generic*mugaz_ref)
     131               RETV_generic = RV_generic/rd_ref-1.
     132               RLvCp_generic = RLVTT_generic/cppd_ref
    137133
    138134               DO k = 1,nlay

  • trunk/LMDZ.GENERIC/libf/phygeneric/thermcell_plume.F90

    r3663 r4146  
    2626      USE tracer_h, ONLY: igcm_h2o_vap
    2727      USE thermcell_mod
    28       USE comcstfi_mod, ONLY: r, cpp, mugaz
    29       USE callkeys_mod, ONLY: water, generic_condensation
     28      USE comcstfi_mod, ONLY: rd_ref, cppd_ref, mugaz_ref
     29      USE callkeys_mod, ONLY: water, generic_condensation, metallicity
    3030      USE generic_cloud_common_h, ONLY: Psat_generic, epsi_generic, RLVTT_generic
    3131
    … …  
    102102      LOGICAL activetmp(ngrid)                        ! If the plume is active (active=true and outgoing mass flux > 0)
    103103
    104       REAL, SAVE :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
    105104      REAL :: RV_generic
    106105      REAL :: RETV_comp, RLvCp_comp !values used for computation (depends if water or generic tracer)
    … …  
    137136      l_start = nlay
    138137
    139       metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
    140 
    141138      ! ALS24 for thermal plume model with generic tracer
    142139      IF (water) THEN
    … …  
    144141         RLvCp_comp = RLvCp
    145142      ELSEIF (generic_condensation .AND. .NOT. water ) THEN
    146          RV_generic = (8.314463*1000.)/(epsi_generic*mugaz)
    147          RETV_comp = RV_generic/r-1.
    148          RLvCp_comp = RLVTT_generic/cpp
     143         RV_generic = (8.314463*1000.)/(epsi_generic*mugaz_ref)
     144         RETV_comp = RV_generic/rd_ref-1.
     145         RLvCp_comp = RLVTT_generic/cppd_ref
    149146      ENDIF
    150147

  • trunk/LMDZ.GENERIC/libf/phygeneric/turbdiff_mod.F90

    r3995 r4146  
    88          ptimestep,pcapcal,                    &   
    99          pplay,pplev,pzlay,pzlev,pz0,                 &
    10           pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf,       &
     10          pu,pv,pt,ph,ppopsk,pq,ptsrf,pemis,pqsurf,       &
    1111          pdtfi,pdqfi,pfluxsrf,            &
    1212          Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
    … …  
    1717      use surfdat_h, only: dryness
    1818      use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
    19       use comcstfi_mod, only: rcp, g, r, cpp
     19      use comcstfi_mod, only: g, rd_ref, cppd_ref
    2020      use callkeys_mod, only: water,tracer,nosurf,kmixmin
    2121      use turb_mod, only : ustar
    … …  
    6161      REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
    6262      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
    63       REAL,INTENT(IN) :: pt(ngrid,nlay),ppopsk(ngrid,nlay)
     63      REAL,INTENT(IN) :: pt(ngrid,nlay),ph(ngrid,nlay),ppopsk(ngrid,nlay)
    6464      REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
    6565      REAL,INTENT(IN) :: pemis(ngrid)
    … …  
    173173         DO ig=1,ngrid
    174174            zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
    175             zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
     175            zExner(ig,ilay)=ppopsk(ig,ilay)
    176176            zovExner(ig,ilay)=1./ppopsk(ig,ilay)
    177177         ENDDO
    178178      ENDDO
    179179
    180       zcst1=4.*g*ptimestep/(R*R)
     180      zcst1=4.*g*ptimestep/(rd_ref*rd_ref)
    181181      DO ilev=2,nlev-1
    182182         DO ig=1,ngrid
    … …  
    186186      ENDDO
    187187      DO ig=1,ngrid
    188          zb0(ig,1)=ptimestep*pplev(ig,1)/(R*ptsrf(ig))
     188         zb0(ig,1)=ptimestep*pplev(ig,1)/(rd_ref*ptsrf(ig))
    189189      ENDDO
    190190      dqsdif_total(:)=0.0
    … …  
    199199            zv(ig,ilev)=pv(ig,ilev)
    200200            zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
    201             zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
     201            zh(ig,ilev)=ph(ig,ilev)!+pdtfi(ig,ilev)*ptimestep*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
    202202        ENDDO
    203203      ENDDO
    … …  
    248248!     ------------------------------------------------------
    249249
    250       call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay,pu,pv,pq,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
    251      
     250      call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay,pplev,pplay, &
     251                   pu,pv,zt,pq,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated
     252                       
    252253!     Adding eddy mixing to mimic 3D general circulation in 1D
    253254!     R. Wordsworth & F. Forget (2010)
    … …  
    273274!               zkv(ig,ilev) = max(kzz_eddy(ilev),zkv(ig,ilev))
    274275!            end do
    275 
    276276            do ig=1,ngrid
    277277               zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
    … …  
    436436
    437437         DO ig=1,ngrid
    438             z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) &
     438            z1(ig)=pcapcal(ig)*ptsrf(ig)+cppd_ref*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) &
    439439                + pfluxsrf(ig)*ptimestep + zdplanck(ig)*ptsrf(ig)
    440             z2(ig) = pcapcal(ig)+zdplanck(ig)+cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))
     440            z2(ig) = pcapcal(ig)+zdplanck(ig)+cppd_ref*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))
    441441            ztsrf(ig) = z1(ig) / z2(ig)
    442442            pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
    … …  
    573573                  zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
    574574
    575                   z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1)   &
     575                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cppd_ref*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1)   &
    576576                      + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
    577577                      + zfluxq(ig,1)*dryness(ig)*RLVTT*((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
    578578
    579                   z2(ig) = pcapcal(ig) + cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
     579                  z2(ig) = pcapcal(ig) + cppd_ref*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
    580580                      + zdplanck(ig)+zfluxq(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*(1.0-zdq(ig,1))
    581581
    … …  
    598598
    599599                      !recompute surface temperature 
    600                       z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1)   &
     600                      z1(ig) = pcapcal(ig)*ptsrf(ig) +cppd_ref*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1)   &
    601601                        + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
    602602                        + RLVTT*dqsdif_total(ig)
    603                       z2(ig) = pcapcal(ig) + cpp*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
     603                      z2(ig) = pcapcal(ig) + cppd_ref*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
    604604                        + zdplanck(ig)
    605605                      ztsrf(ig) = z1(ig) / z2(ig)
    … …  
    720720     
    721721      DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
    722          sensibFlux(ig)=cpp*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
     722         sensibFlux(ig)=cppd_ref*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
    723723      ENDDO
    724724

  • trunk/LMDZ.GENERIC/libf/phygeneric/vdif_kc.F

    r3662 r4146  
    1       SUBROUTINE vdif_kc(ngrid,nlay,nq,dt,g,zlev,zlay,u,v,
    2      &                       zq,teta,cd,q2,km,kn)
     1      SUBROUTINE vdif_kc(ngrid,nlay,nq,dt,g,zlev,zlay,pplev,pplay,u,v,
     2     &                       zt,zq,teta,cd,q2,km,kn)
    33      use generic_cloud_common_h, only: epsi_generic
    44      use generic_tracer_index_mod, only: generic_tracer_index
    55      use callkeys_mod, only: generic_condensation,
    6      &                        virtual_theta_correction
     6     &                        virtual_theta_correction,thermo_phy
     7      use thermo_mod, only: rcp
     8      use tracer_h
    79      IMPLICIT NONE
    810c.......................................................................
    … …  
    3537      REAL,INTENT(IN) :: zlev(ngrid,nlay+1)
    3638      REAL,INTENT(IN) :: zlay(ngrid,nlay)
     39      REAL,INTENT(IN) :: pplev(ngrid,nlay+1)
     40      REAL,INTENT(IN) :: pplay(ngrid,nlay)
    3741      REAL,INTENT(IN) :: u(ngrid,nlay)
    3842      REAL,INTENT(IN) :: v(ngrid,nlay)
     43      REAL,INTENT(IN) :: zt(ngrid,nlay)
    3944      REAL,INTENT(IN) :: zq(ngrid,nlay,nq)
    4045      REAL,INTENT(IN) :: teta(ngrid,nlay)
    … …  
    5964      REAL unsdzdec(ngrid,nlay+1)
    6065      REAL q(ngrid,nlay+1)
     66      REAL ztv(ngrid,nlay)
    6167      REAL tetav(ngrid,nlay)
    6268      integer iq,igcm_generic_vap, igcm_generic_ice
    … …  
    263269c
    264270c.......................................................................
     271c  Brunt-Vaisala frequency
     272c.......................................................................
     273c
     274
     275      SELECT CASE (TRIM(thermo_phy))
     276     
     277      CASE ('thermo_uni_ideal')
     278
     279
    265280c  Virtual theta correction
    266 c.......................................................................
    267 c
    268281
    269282      if((generic_condensation) .and. (virtual_theta_correction)) THEN
    … …  
    283296      endif
    284297
    285 c
    286 c-----------------------------------------------------------------------
    287298      DO ilev=2,nlev-1
    288299       DO igrid=1,ngrid
    289 c-----------------------------------------------------------------------
    290 c
     300
    291301      if((generic_condensation) .and. (virtual_theta_correction)) THEN
    292302        n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
    … …  
    317327        m(igrid,ilev)=sqrt(m2(igrid,ilev))
    318328        mpre(igrid,ilev)=m(igrid,ilev)
    319 c
    320 c-----------------------------------------------------------------------
     329
    321330       ENDDO
    322331      ENDDO
    323 c-----------------------------------------------------------------------
    324 c
     332     
     333      END SELECT
     334c-----------------------------------------------------------------------
     335c
     336
     337
    325338      DO igrid=1,ngrid
    326339        m2(igrid,nlev)=m2(igrid,nlev-1)

  • trunk/LMDZ.GENERIC/libf/phygeneric/vdifc_mod.F

    r3236 r4146  
    88     &     ptimestep,pcapcal,lecrit,                       
    99     &     pplay,pplev,pzlay,pzlev,pz0,
    10      &     pu,pv,ph,pq,ptsrf,pemis,pqsurf,
     10     &     pu,pv,pt,ph,pq,ptsrf,pemis,pqsurf,
    1111     &     pdhfi,pdqfi,pfluxsrf,
    1212     &     pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
    … …  
    1818      use surfdat_h, only: dryness
    1919      use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
    20       use comcstfi_mod, only: g, r, cpp, rcp
     20      use comcstfi_mod, only: g, rd_ref, cppd_ref, rcp_ref
    2121      use callkeys_mod, only: water,tracer,nosurf
    2222
    … …  
    5353      REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    5454      REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
    55       REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
     55      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
     56      REAL,INTENT(IN) :: pt(ngrid,nlay),ph(ngrid,nlay)
    5657      REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid)
    5758      REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
    … …  
    179180      ENDDO
    180181
    181       zcst1=4.*g*ptimestep/(R*R)
     182      zcst1=4.*g*ptimestep/(rd_ref*rd_ref)
    182183      DO ilev=2,nlev-1
    183184         DO ig=1,ngrid
    184185            zb0(ig,ilev)=pplev(ig,ilev)*
    185      s           (pplev(ig,1)/pplev(ig,ilev))**rcp /
     186     s           (pplev(ig,1)/pplev(ig,ilev))**rcp_ref /
    186187     s           (ph(ig,ilev-1)+ph(ig,ilev))
    187188            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
    … …  
    190191      ENDDO
    191192      DO ig=1,ngrid
    192          zb0(ig,1)=ptimestep*pplev(ig,1)/(R*ptsrf(ig))
     193         zb0(ig,1)=ptimestep*pplev(ig,1)/(rd_ref*ptsrf(ig))
    193194      ENDDO
    194195
    … …  
    246247!     ------------------------------------------------------
    247248
    248       call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
    249      &     ,pu,pv,pq,ph,zcdv_true
     249      call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay,
     250     &     pplev,pplay,pu,pv,pt,pq,ph,zcdv_true
    250251     &     ,pq2,zkv,zkh)
    251252
    … …  
    405406         DO ig=1,ngrid
    406407
    407             z1(ig) = pcapcal(ig)*ptsrf(ig) + cpp*zb(ig,1)*zc(ig,1)
    408      &           + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep
    409             z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1))
    410      &           +zdplanck(ig)
    411             ztsrf2(ig) = z1(ig) / z2(ig)
    412             pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
    413             zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
     408           z1(ig)=pcapcal(ig)*ptsrf(ig)+cppd_ref*zb(ig,1)*zc(ig,1)
     409     &          + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep
     410           z2(ig) = pcapcal(ig) + cppd_ref*zb(ig,1)*(1.-zd(ig,1))
     411     &          +zdplanck(ig)
     412           ztsrf2(ig) = z1(ig) / z2(ig)
     413           pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
     414           zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
    414415         ENDDO
    415416
    … …  
    552553! calculation of h0 and h1
    553554               do ig=1,ngrid
    554                   zdq0(ig) = dqsat(ig)
    555                   zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
    556 
    557                   z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zb(ig,1)*zc(ig,1)
    558      &                 + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep
    559      &                 +  zb(ig,1)*dryness(ig)*RLVTT*
    560      &                 ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
    561 
    562                   z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1))
    563      &                 +zdplanck(ig)
    564      &                 +zb(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*
    565      &                 (1.0-zdq(ig,1))
    566 
    567                   ztsrf2(ig) = z1(ig) / z2(ig)
    568                   pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
    569                   zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
     555                 zdq0(ig) = dqsat(ig)
     556                 zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
     557
     558                 z1(ig)=pcapcal(ig)*ptsrf(ig)+cppd_ref*zb(ig,1)*zc(ig,1)
     559     &                + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep
     560     &                +  zb(ig,1)*dryness(ig)*RLVTT*
     561     &                ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
     562
     563                 z2(ig) = pcapcal(ig) + cppd_ref*zb(ig,1)*(1.-zd(ig,1))
     564     &                +zdplanck(ig)
     565     &                +zb(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*
     566     &                (1.0-zdq(ig,1))
     567
     568                 ztsrf2(ig) = z1(ig) / z2(ig)
     569                 pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
     570                 zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
    570571               enddo
    571572
    … …  
    664665
    665666      DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
    666          sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zh(ig,1)-ztsrf2(ig))
     667        sensibFlux(ig)=cppd_ref*zb(ig,1)/ptimestep*(zh(ig,1)-ztsrf2(ig))
    667668      ENDDO     
    668669

  • trunk/LMDZ.GENERIC/libf/phygeneric/watercommon_h.F90

    r3663 r4146  
    2828      subroutine su_watercycle
    2929
    30       use comcstfi_mod, only: r, cpp, mugaz
     30      use comcstfi_mod, only: rd_ref, cppd_ref, mugaz_ref
    3131      implicit none
    3232
    … …  
    4545!==================================================================
    4646
    47       epsi   = mH2O / mugaz
    48       RCPD   = cpp 
     47      epsi   = mH2O / mugaz_ref
     48      RCPD   = cppd_ref
    4949
    5050      !RW = 1000.*R/mH2O
    … …  
    5656     
    5757! AB : initializations added for the thermal plume model
    58       RETV = RW / r - 1.
     58      RETV = RW / rd_ref - 1.
    5959      RLvCp = RLVTT / RCPD
    6060     
Note: See TracChangeset for help on using the changeset viewer.