Changeset 4982 for LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

Timestamp:

Jun 15, 2024, 5:17:08 PM (11 days ago)

Author:

crisi

Message:

suppress isotope_params.def + update physiq_mod + proof of concept of 3rd dimension with reevap routine

File:

• LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90 (modified) (11 diffs)

LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

@@ -20 +20 @@
 
    !--- Variables not depending on isotopes
-   REAL,    SAVE :: pxtmelt, pxtice, pxtmin, pxtmax
-!$OMP THREADPRIVATE(pxtmelt, pxtice, pxtmin, pxtmax)
-   REAL,    SAVE :: tdifexp, tv0cin, thumxt1
-!$OMP THREADPRIVATE(tdifexp, tv0cin, thumxt1)
+   REAL,    SAVE :: thumxt1
+   !$OMP THREADPRIVATE(thumxt1)
    INTEGER, SAVE :: ntot 
 !$OMP THREADPRIVATE(ntot)
@@ -30 +28 @@
    REAL,    SAVE :: P_veg
 !$OMP THREADPRIVATE(P_veg)
-   REAL,    SAVE :: musi, lambda_sursat
-!$OMP THREADPRIVATE(musi, lambda_sursat)
-   REAL,    SAVE :: Kd
-!$OMP THREADPRIVATE(Kd)
-   REAL,    SAVE :: rh_cste_surf_cond, T_cste_surf_cond
-!$OMP THREADPRIVATE(rh_cste_surf_cond, T_cste_surf_cond)
+   REAL,    SAVE :: lambda_sursat
+!$OMP THREADPRIVATE(lambda_sursat)
    LOGICAL, SAVE :: bidouille_anti_divergence    ! T: regularly, xteau <- q to avoid slow drifts
 !$OMP THREADPRIVATE(bidouille_anti_divergence)
@@ -54 +48 @@
    REAL,    SAVE :: sstlatcrit, dsstlatcrit
 !$OMP THREADPRIVATE(sstlatcrit, dsstlatcrit)
-   REAL,    SAVE :: deltaO18_oce
-!$OMP THREADPRIVATE(deltaO18_oce)
    INTEGER, SAVE :: albedo_prescrit              ! 0: default ; 1: constant albedo
 !$OMP THREADPRIVATE(albedo_prescrit)
@@ -88 +80 @@
    REAL,    SAVE :: fac_modif_evaoce
 !$OMP THREADPRIVATE(fac_modif_evaoce)
+   REAL,    SAVE :: deltaO18_oce
+!$OMP THREADPRIVATE(deltaO18_oce)
    INTEGER, SAVE :: ok_bidouille_wake
 !$OMP THREADPRIVATE(ok_bidouille_wake)
@@ -106 +100 @@
                     alpha_liq_sol, Rdefault, Rmethox
 !$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox)
-   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
-!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
+!   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
+!!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
+
+   !--- H2[18]O reference
+   REAL, PARAMETER :: fac_enrichoce18=0.0005
+   REAL, PARAMETER :: alpha_liq_sol_O18=1.00291
+   REAL, PARAMETER :: talph1_O18=1137.
+   REAL, PARAMETER :: talph2_O18=-0.4156
+   REAL, PARAMETER :: talph3_O18=-2.0667E-3
+   REAL, PARAMETER :: talps1_O18=11.839
+   REAL, PARAMETER :: talps2_O18=-0.028244
+   REAL, PARAMETER :: tdifrel_O18=1./0.9723
+   REAL, PARAMETER :: tkcin0_O18=0.006
+   REAL, PARAMETER :: tkcin1_O18=0.000285
+   REAL, PARAMETER :: tkcin2_O18=0.00082
+   REAL, PARAMETER :: fac_coeff_eq17_liq=0.529
+   REAL, PARAMETER :: fac_coeff_eq17_ice=0.529
+
+   !---- Parameters that do not depend on the nature of water isotopes:
+   REAL, PARAMETER :: pxtmelt = 273.15 ! temperature at which ice formation starts
+   REAL, PARAMETER :: pxtice  = 273.15-10.0 ! -- temperature at which all condensate is ice:
+   REAL, PARAMETER :: pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
+   REAL, PARAMETER :: pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
+   REAL, PARAMETER :: tdifexp = 0.58 ! -- a constant for alpha_eff for equilibrium below cloud base:
+   REAL, PARAMETER :: tv0cin  = 7.0 ! wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979)
+   REAL, PARAMETER :: musi=1.0  ! facteurs lambda et mu dans Si=musi-lambda*T
+   REAL, PARAMETER :: Kd=2.5e-9 ! m2/s ! diffusion dans le sol
+   REAL, PARAMETER :: rh_cste_surf_cond = 0.6 ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
+   REAL, PARAMETER :: T_cste_surf_cond = 288.0
+
 
    !--- Negligible lower thresholds: no need to check for absurd values under these lower limits
@@ -140 +162 @@
    INTEGER :: ixt
 
-   !--- H2[18]O reference
-   REAL :: fac_enrichoce18, alpha_liq_sol_O18, &
-           talph1_O18, talph2_O18, talph3_O18, talps1_O18, talps2_O18, &
-           tkcin0_O18, tkcin1_O18, tkcin2_O18, tdifrel_O18  
-
+ 
    !--- For H2[17]O
    REAL    :: fac_kcin, pente_MWL
@@ -152 +170 @@
    LOGICAL, PARAMETER ::   ok_nocinsat = .FALSE. ! if T: no sursaturation effect for ice
    LOGICAL, PARAMETER :: Rdefault_smow = .FALSE. ! if T: Rdefault=smow; if F: nul
+   LOGICAL, PARAMETER :: tnat1 = .TRUE. ! If T: all tnats are 1.
 
    !--- For [3]H
@@ -157 +176 @@
 
    CHARACTER(LEN=maxlen) :: modname, sxt
-   REAL, ALLOCATABLE :: tmp(:)
 
    modname = 'iso_init'
@@ -214 +232 @@
       CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
    END IF
-   deltaO18_oce=0.0
    CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
    CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
@@ -249 +266 @@
    CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
 
-   !--------------------------------------------------------------
-   ! Parameters that do not depend on the nature of water isotopes:
-   !--------------------------------------------------------------
-   ! -- temperature at which ice condensate starts to form (valeur ECHAM?):
-   pxtmelt = 273.15
-
-   ! -- temperature at which all condensate is ice:
-   pxtice  = 273.15-10.0
-
-   !- -- test PHASE
-!   pxtmelt = 273.15 - 10.0
-!   pxtice  = 273.15 - 30.0
-
-   ! -- minimum temperature to calculate fractionation coeff
-   pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
-   pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
-   !    Remarque: les coeffs ont ete mesures seulement jusq'à -40!
-
-   ! -- a constant for alpha_eff for equilibrium below cloud base:
-   tdifexp = 0.58
-   tv0cin  = 7.0
-
-   ! facteurs lambda et mu dans Si=musi-lambda*T
-   musi=1.0
-   if (ok_nocinsat) lambda_sursat = 0.0          ! no sursaturation effect
-
-   ! diffusion dans le sol
-   Kd=2.5e-9 ! m2/s   
-
-   ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
-   rh_cste_surf_cond = 0.6
-    T_cste_surf_cond = 288.0
+   ! Ocean composition
+   CALL get_in('deltaO18_oce',  deltaO18_oce, 0.0)
    
    CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
 
    !--------------------------------------------------------------
-   ! Parameters that depend on the nature of water isotopes:
+   ! Isotope fractionation factors and a few isotopic constants
    !--------------------------------------------------------------
-   IF(getKey('tnat',    tnat,    isoName)) CALL abort_physic(modname, 'can''t get tnat',    1)
-   IF(getKey('toce',    toce,    isoName)) CALL abort_physic(modname, 'can''t get toce',    1)
-   IF(getKey('tcorr',   tcorr,   isoName)) CALL abort_physic(modname, 'can''t get tcorr',   1)
-   IF(getKey('talph1',  talph1,  isoName)) CALL abort_physic(modname, 'can''t get talph1',  1)
-   IF(getKey('talph2',  talph2,  isoName)) CALL abort_physic(modname, 'can''t get talph2',  1)
-   IF(getKey('talph3',  talph3,  isoName)) CALL abort_physic(modname, 'can''t get talph3',  1)
-   IF(getKey('talps1',  talps1,  isoName)) CALL abort_physic(modname, 'can''t get talps1',  1)
-   IF(getKey('talps2',  talps2,  isoName)) CALL abort_physic(modname, 'can''t get talps2',  1)
-   IF(getKey('tkcin0',  tkcin0,  isoName)) CALL abort_physic(modname, 'can''t get tkcin0',  1)
-   IF(getKey('tkcin1',  tkcin1,  isoName)) CALL abort_physic(modname, 'can''t get tkcin1',  1)
-   IF(getKey('tkcin2',  tkcin2,  isoName)) CALL abort_physic(modname, 'can''t get tkcin2',  1)
-   IF(getKey('tdifrel', tdifrel, isoName)) CALL abort_physic(modname, 'can''t get tdifrel', 1)
-   IF(getKey('alpha_liq_sol', alpha_liq_sol, isoName)) CALL abort_physic(modname, 'can''t get alpha_liq_sol',  1)
-   IF(getKey('Rdefault',Rdefault,isoName)) CALL abort_physic(modname, 'can''t get Rdefault',1)
-   IF(getKey('Rmethox', Rmethox, isoName)) CALL abort_physic(modname, 'can''t get Rmethox', 1)
+   ALLOCATE(tkcin0(niso))
+   ALLOCATE(tkcin1(niso))
+   ALLOCATE(tkcin2(niso))
+   ALLOCATE(tnat(niso))
+   ALLOCATE(tdifrel(niso))
+   ALLOCATE(toce(niso))
+   ALLOCATE(tcorr(niso))
+   ALLOCATE(talph1(niso))
+   ALLOCATE(talph2(niso))
+   ALLOCATE(talph3(niso))
+   ALLOCATE(talps1(niso))
+   ALLOCATE(talps2(niso))
+   ALLOCATE(alpha_liq_sol(niso))
+   ALLOCATE(Rdefault(niso))
+   ALLOCATE(Rmethox(niso))
+
+   do ixt=1,niso
+     if (ixt.eq.iso_HTO) then  ! Tritium
+       tkcin0(ixt) = 0.01056
+       tkcin1(ixt) = 0.0005016
+       tkcin2(ixt) = 0.0014432
+       if (tnat1) then 
+               tnat(ixt)=1 
+       else 
+               tnat(ixt)=0.
+       endif
+       toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
+       tcorr(ixt)=1.
+       tdifrel(ixt)=1./0.968
+       talph1(ixt)=46480.
+       talph2(ixt)=-103.87
+       talph3(ixt)=0.
+       talps1(ixt)=46480.
+       talps2(ixt)=-103.87
+       alpha_liq_sol(ixt)=1. 
+       Rmethox(ixt)=0.0
+     endif
+     if (ixt.eq.iso_O17) then  ! O17
+       pente_MWL=0.528 
+       tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG ! tdifrel(ixt)=1./0.985452 ! donné par Amaelle
+       fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145 ! donné par Amaelle
+       tkcin0(ixt) = tkcin0_O18*fac_kcin
+       tkcin1(ixt) = tkcin1_O18*fac_kcin
+       tkcin2(ixt) = tkcin2_O18*fac_kcin
+       if (tnat1) then 
+               tnat(ixt)=1 
+       else 
+               tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
+       endif
+       toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
+       tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle            
+       talph1(ixt)=talph1_O18
+       talph2(ixt)=talph2_O18
+       talph3(ixt)=talph3_O18
+       talps1(ixt)=talps1_O18
+       talps2(ixt)=talps2_O18     
+       alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq 
+       Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
+       Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
+     endif
+     if (ixt.eq.iso_O18) then  ! Oxygene18
+       tkcin0(ixt) = tkcin0_O18
+       tkcin1(ixt) = tkcin1_O18
+       tkcin2(ixt) = tkcin2_O18
+       if (tnat1) then
+               tnat(ixt)=1 
+       else 
+               tnat(ixt)=2005.2E-6
+       endif
+       toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
+       tcorr(ixt)=1.0+fac_enrichoce18
+       tdifrel(ixt)=tdifrel_O18
+       talph1(ixt)=talph1_O18
+       talph2(ixt)=talph2_O18
+       talph3(ixt)=talph3_O18
+       talps1(ixt)=talps1_O18
+       talps2(ixt)=talps2_O18
+       alpha_liq_sol(ixt)=alpha_liq_sol_O18   
+       Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
+       Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006  
+     endif
+     if (ixt.eq.iso_HDO) then ! Deuterium
+       pente_MWL=8.0
+       tdifrel(ixt)=1./0.9755 !          fac_kcin=0.88
+       fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
+       tkcin0(ixt) = tkcin0_O18*fac_kcin
+       tkcin1(ixt) = tkcin1_O18*fac_kcin
+       tkcin2(ixt) = tkcin2_O18*fac_kcin
+       if (tnat1) then 
+               tnat(ixt)=1
+       else
+               tnat(ixt)=155.76E-6
+       endif
+       toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
+       tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL          
+       talph1(ixt)=24844.
+       talph2(ixt)=-76.248
+       talph3(ixt)=52.612E-3
+       talps1(ixt)=16288.
+       talps2(ixt)=-0.0934
+       !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
+       alpha_liq_sol(ixt)=1.0212 
+       ! valeur de Lehmann & Siegenthaler, 1991, Journal of
+       ! Glaciology, vol 37, p 23
+       Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
+       Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
+     endif
+     if (ixt.eq.iso_eau) then ! Oxygene16
+       tkcin0(ixt) = 0.0
+       tkcin1(ixt) = 0.0
+       tkcin2(ixt) = 0.0
+       tnat(ixt)=1.
+       toce(ixt)=tnat(ixt)
+       tcorr(ixt)=1.0
+       tdifrel(ixt)=1.
+       talph1(ixt)=0.
+       talph2(ixt)=0.
+       talph3(ixt)=0.
+       talps1(ixt)=0.
+       talph3(ixt)=0.
+       alpha_liq_sol(ixt)=1.
+       Rdefault(ixt)=tnat(ixt)*1.0
+       Rmethox(ixt)=1.0
+     endif
+   enddo ! ixt=1,niso
 
    IF(.NOT.Rdefault_smow) then
@@ -308 +405 @@
    ENDIF
    write(*,*) 'Rdefault=',Rdefault
+   write(*,*) 'toce=',toce
 
    !--- Sensitivity test: no kinetic effect in sfc evaporation