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

Timestamp:

May 19, 2022, 10:47:23 AM (2 years ago)

Author:

dcugnet

Message:

• include new type "tracer.def" tracers description files in DefLists/?
• Remove the old code (lOldCode=.TRUE.)
• Remove the hard-coded part in isotopes_mod that was ensuring the convergence => this version will not converge with previous isotopic physics, but the

difference, only due to round off errors on "tdifrel", is acceptable.

File:

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

LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

@@ -11 +11 @@
   !--- Contains all isotopic variables + their initialization
   !--- Isotopes-specific routines are in isotopes_routines_mod to avoid circular dependencies with isotopes_verif_mod.
-
-  LOGICAL, PARAMETER :: lOldCode=.FALSE.
 
    !--- Isotopes indices (in [1,niso] ; non-existing => 0 index)
@@ -165 +163 @@
    CALL msg('219: entree', modname)
 
-   !--- Memory allocations
-   IF(lOldCode) THEN
-   ALLOCATE(talph1(niso), tkcin0(niso),  talps1(niso),  tnat(niso))
-   ALLOCATE(talph2(niso), tkcin1(niso),  talps2(niso),  toce(niso))
-   ALLOCATE(talph3(niso), tkcin2(niso), tdifrel(niso), tcorr(niso))
-   ALLOCATE(alpha_liq_sol(niso),   Rdefault(niso),   Rmethox(niso))
-   END IF
-
-
    !--------------------------------------------------------------
    ! General:
@@ -231 +220 @@
       CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
    END IF
-   IF(lOldCode) &
-   CALL get_in('deltaO18_oce',        deltaO18_oce,   0.0)
    deltaO18_oce=0.0
    CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
@@ -306 +293 @@
    ! Parameters that depend on the nature of water isotopes:
    !--------------------------------------------------------------
-
-   !===========================================================================================================================
-   IF(lOldCode) THEN
-   !===========================================================================================================================
-
-   ! Local constants
-   fac_enrichoce18 = 0.0005            ! Then: tcorO18 = 1 + fac_enrichoce18
-                                       !       tcorD   = 1 + fac_enrichoce18*8
-                                       !       tcorO17 = 1 + fac_enrichoce18*0.528
-   alpha_liq_sol_O18 = 1.00291         ! From Lehmann & Siegenthaler, 1991, 
-                                       ! Journal of Glaciology, vol 37, p 23
-   talph1_O18 = 1137.  ; talph2_O18 = -0.4156   ; talph3_O18 = -2.0667E-3
-   talps1_O18 = 11.839 ; talps2_O18 = -0.028244 
-   tkcin0_O18 = 0.006  ; tkcin1_O18 =  0.000285 ; tkcin2_O18 =  0.00082
-   tdifrel_O18 = 1./0.9723
-
-   ! ln(alphaeq) ratio between O18 and O17
-   fac_coeff_eq17_liq = 0.529          ! From Amaelle
-  !fac_coeff_eq17_ice = 0.528          ! slope MWL 
-   fac_coeff_eq17_ice = 0.529
-
-   !--- Kinetic factor for surface evaporation:
-   ! (cf: kcin = tkcin0                  if |V|<tv0cin
-   !      kcin = tkcin1*|Vsurf| + tkcin2 if |V|>tv0cin )
-   ! (Rq: formula discontinuous for |V|=tv0cin... )        
-
-   DO ixt = 1, niso
-      sxt=int2str(ixt)
-      CALL msg('80: ixt='//TRIM(int2str(ixt)),modname)
-
-      Rdefault(ixt) = 0.0
-      IF(ixt == iso_HTO) THEN          !=== H[3]HO
-         tdifrel(ixt) = 1./0.968
-         tkcin0(ixt) = 0.01056
-         tkcin1(ixt) = 0.0005016
-         tkcin2(ixt) = 0.0014432
-         tnat  (ixt) = 0.
-         toce  (ixt) = 4.0E-19         ! Ratio T/H = 0.2 TU, Dreisigacker and Roether 1978
-        !toce  (ixt) = 2.2222E-8       ! Corrected by Alex Cauquoin
-        !toce  (ixt) = 1.0E-18         ! Ratio 3H/1H ocean
-         tcorr (ixt) = 1.
-         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
-      ELSE IF(ixt == iso_O17) THEN     !=== H2[17]O
-         tdifrel(ixt)=1./0.98555       ! Used in 1D and in LdG's model
-        !tdifrel(ixt)=1./0.985452      ! From Amaelle 
-        !fac_kcin=0.5145               ! From Amaelle
-         fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
-         tkcin0(ixt) = tkcin0_O18*fac_kcin
-         tkcin1(ixt) = tkcin1_O18*fac_kcin
-         tkcin2(ixt) = tkcin2_O18*fac_kcin
-         tnat  (ixt) = 0.004/100.      ! O17 = 0.004% of oxygen
-         pente_MWL=0.528
-         toce  (ixt) = tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
-         tcorr (ixt) = 1.0+fac_enrichoce18*pente_MWL  ! From 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 
-         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*(-3.15/1000.0+1.0)
-         Rmethox(ixt) = (230./1000.+1.)*tnat(ixt)     ! Zahn et al 2006
-      ELSE IF(ixt == iso_O18) THEN     !=== H2[18]O
-         tdifrel(ixt) = tdifrel_O18
-         tkcin0(ixt) = tkcin0_O18
-         tkcin1(ixt) = tkcin1_O18
-         tkcin2(ixt) = tkcin2_O18
-         tnat  (ixt) = 2005.2E-6
-         toce  (ixt) = tnat(ixt)*(1.0+deltaO18_oce/1000.0)
-         tcorr (ixt) = 1.0+fac_enrichoce18
-         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   
-         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*(-6.0/1000.0+1.0)
-         Rmethox(ixt) = (130./1000.+1.)*tnat(ixt) ! Zahn et al 2006   
-      ELSE IF(ixt == iso_HDO) THEN     !=== H[2]HO
-         tdifrel(ixt) = 1./0.9755
-        !fac_kcin=0.88
-         fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
-         tkcin0(ixt) = tkcin0_O18*fac_kcin
-         tkcin1(ixt) = tkcin1_O18*fac_kcin
-         tkcin2(ixt) = tkcin2_O18*fac_kcin
-         tnat  (ixt) = 155.76E-6
-         pente_MWL = 8.0
-         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
-        !alpha_liq_sol(ixt)=1.0192 ZX  ! From Weston, Ralph, 1955
-         alpha_liq_sol(ixt)=1.0212     ! From Lehmann & Siegenthaler, 1991, 
-                                       ! Journal of Glaciology, vol 37, p 23
-         IF(Rdefault_smow) 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
-      ELSE IF(ixt  == iso_eau) THEN    !=== H2O[16]
-         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. ; talps2(ixt) = 0.
-         alpha_liq_sol(ixt)=1.
-         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*1.0
-         Rmethox(ixt) = 1.0
-      END IF
-   END DO
-   !===========================================================================================================================
-   ELSE
-   !===========================================================================================================================
-
    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)
@@ -430 +305 @@
    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)
-   DO ixt = 1, niso
-      IF     (ixt == iso_HTO) THEN; tdifrel(ixt) = 1./0.968
-      ELSE IF(ixt == iso_HDO) THEN; tdifrel(ixt) = 1./0.9755
-      ELSE IF(ixt == iso_O17) THEN; tdifrel(ixt) = 1./0.98555
-      ELSE IF(ixt == iso_O18) THEN; tdifrel(ixt) = 1./0.9723
-      ELSE IF(ixt == iso_eau) THEN; tdifrel(ixt) = 1.; END IF
-   END DO
    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)
    IF(.NOT.Rdefault_smow) Rdefault(:) = 0.0
-
-   !===========================================================================================================================
-   END IF
-   !===========================================================================================================================
 
    !--- Sensitivity test: no kinetic effect in sfc evaporation