Ignore:
Timestamp:
May 14, 2022, 8:13:22 PM (2 years ago)
Author:
dcugnet
Message:
  • remove striso (use isoName instead)
  • few fixes for the lOldCode=.FALSE. code
  • add the « isotopes_params.def » file, used in the lOldCode=.FALSE. part of the isotopes_mod module.


File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

    r4143 r4149  
    33
    44MODULE isotopes_mod
    5    USE strings_mod, ONLY: msg, real2str, int2str, bool2str, maxlen, strIdx, strStack
     5   USE strings_mod,  ONLY: msg, real2str, int2str, bool2str, maxlen, strIdx, strStack
     6   USE infotrac_phy, ONLY: isoName
    67   IMPLICIT NONE
    78   INTERFACE get_in; MODULE PROCEDURE getinp_s, getinp_i, getinp_r, getinp_l;  END INTERFACE get_in
    … …  
    1011  !--- Contains all isotopic variables + their initialization
    1112  !--- Isotopes-specific routines are in isotopes_routines_mod to avoid circular dependencies with isotopes_verif_mod.
     13
     14  LOGICAL, PARAMETER :: lOldCode=.TRUE.
    1215
    1316   !--- Isotopes indices (in [1,niso] ; non-existing => 0 index)
    … …  
    105108                    alpha_liq_sol, Rdefault, Rmethox
    106109!$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox)
    107 character*3, ALLOCATABLE, DIMENSION(:), save :: striso
    108 !$OMP THREADPRIVATE(striso)
    109110   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
    110111!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
    … …  
    136137SUBROUTINE iso_init()
    137138   USE ioipsl_getin_p_mod, ONLY: getin_p
    138    USE infotrac_phy,       ONLY: ntiso, niso, isoName
     139   USE infotrac_phy,       ONLY: ntiso, niso, getKey
     140    USE strings_mod,       ONLY: maxlen
    139141   IMPLICIT NONE
    140142
    … …  
    149151   !--- For H2[17]O
    150152   REAL    :: fac_kcin, pente_MWL
    151    INTEGER :: ierr
    152153     
    153154   !--- Sensitivity tests
    … …  
    160161
    161162   CHARACTER(LEN=maxlen) :: modname, sxt
     163   REAL, ALLOCATABLE :: tmp(:)
    162164
    163165   modname = 'iso_init'
    … …  
    165167
    166168   !--- Memory allocations
     169   IF(lOldCode) THEN
    167170   ALLOCATE(talph1(niso), tkcin0(niso),  talps1(niso),  tnat(niso))
    168171   ALLOCATE(talph2(niso), tkcin1(niso),  talps2(niso),  toce(niso))
    169172   ALLOCATE(talph3(niso), tkcin2(niso), tdifrel(niso), tcorr(niso))
    170173   ALLOCATE(alpha_liq_sol(niso),   Rdefault(niso),   Rmethox(niso))
    171    ALLOCATE(striso(niso))
     174   END IF
    172175
    173176
    … …  
    184187
    185188   !--- Type of water isotopes:
    186    iso_eau = strIdx(isoName, 'H2[16]O'); CALL msg('59: iso_eau='//int2str(iso_eau), modname)
    187    iso_O17 = strIdx(isoName, 'H2[17]O'); CALL msg('iso_HDO='//int2str(iso_HDO), modname)
     189   iso_eau = strIdx(isoName, 'H2[16]O'); CALL msg('iso_eau='//int2str(iso_eau), modname)
     190   iso_HDO = strIdx(isoName, 'H[2]HO'); CALL msg('iso_HDO='//int2str(iso_HDO), modname)
    188191   iso_O18 = strIdx(isoName, 'H2[18]O'); CALL msg('iso_O18='//int2str(iso_O18), modname)
    189    iso_HDO = strIdx(isoName, 'H[2]HO'); CALL msg('iso_O17='//int2str(iso_O17), modname)
     192   iso_O17 = strIdx(isoName, 'H2[17]O'); CALL msg('iso_O17='//int2str(iso_O17), modname)
    190193   iso_HTO = strIdx(isoName, 'H[3]HO');  CALL msg('iso_HTO='//int2str(iso_HTO), modname)
    191194
    … …  
    229232      CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
    230233   END IF
     234   IF(lOldCode) &
    231235   CALL get_in('deltaO18_oce',        deltaO18_oce,   0.0)
    232 
     236   deltaO18_oce=0.0
    233237   CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
    234238   CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
    … …  
    298302    T_cste_surf_cond = 288.0
    299303   
     304   CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
     305
    300306   !--------------------------------------------------------------
    301307   ! Parameters that depend on the nature of water isotopes:
    302308   !--------------------------------------------------------------
     309
     310   !===========================================================================================================================
     311   IF(lOldCode) THEN
     312   !===========================================================================================================================
     313
    303314   ! Local constants
    304315   fac_enrichoce18 = 0.0005            ! Then: tcorO18 = 1 + fac_enrichoce18
    … …  
    317328   fac_coeff_eq17_ice = 0.529
    318329
    319    CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
    320 
    321330   !--- Kinetic factor for surface evaporation:
    322331   ! (cf: kcin = tkcin0                  if |V|<tv0cin
    … …  
    326335   DO ixt = 1, niso
    327336      sxt=int2str(ixt)
    328       WRITE(*,*) 'iso_init 80: ixt=',ixt
     337      CALL msg('80: ixt='//TRIM(int2str(ixt)),modname)
    329338
    330339      Rdefault(ixt) = 0.0
    … …  
    343352         alpha_liq_sol(ixt) = 1.
    344353         Rmethox(ixt) = 0.0
    345          striso (ixt) = 'HTO'
    346354      ELSE IF(ixt == iso_O17) THEN     !=== H2[17]O
    347355         tdifrel(ixt)=1./0.98555       ! Used in 1D and in LdG's model
    … …  
    361369         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*(-3.15/1000.0+1.0)
    362370         Rmethox(ixt) = (230./1000.+1.)*tnat(ixt)     ! Zahn et al 2006
    363          striso (ixt) = 'O17'
    364371      ELSE IF(ixt == iso_O18) THEN     !=== H2[18]O
    365372         tdifrel(ixt) = tdifrel_O18
    … …  
    375382         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*(-6.0/1000.0+1.0)
    376383         Rmethox(ixt) = (130./1000.+1.)*tnat(ixt) ! Zahn et al 2006   
    377          striso (ixt) = 'O18'
    378          CALL msg('519: ixt, striso(ixt) = '//TRIM(sxt)//', '//TRIM(striso(ixt)), modname)
    379384      ELSE IF(ixt == iso_HDO) THEN     !=== H[2]HO
    380385         tdifrel(ixt) = 1./0.9755
    … …  
    395400         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
    396401         Rmethox(ixt) = tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
    397          striso (ixt) = 'HDO'
    398          CALL msg('548: ixt,striso(ixt) = '//TRIM(sxt)//', '//striso(ixt), modname)
    399402      ELSE IF(ixt  == iso_eau) THEN    !=== H2O[16]
    400403         tkcin0(ixt) = 0.0
    … …  
    406409         tdifrel(ixt) = 1.
    407410         talph1(ixt) = 0. ; talph2(ixt) = 0. ; talph3(ixt) = 0.
    408          talps1(ixt) = 0. ; talph3(ixt) = 0.
     411         talps1(ixt) = 0. ; talps2(ixt) = 0.
    409412         alpha_liq_sol(ixt)=1.
    410413         IF(Rdefault_smow) Rdefault(ixt) = tnat(ixt)*1.0
    411414         Rmethox(ixt) = 1.0
    412          striso(ixt) = 'eau'
    413415      END IF
    414416   END DO
     417   !===========================================================================================================================
     418   ELSE
     419   !===========================================================================================================================
     420
     421   IF(getKey('tnat',    tnat,    isoName)) CALL abort_physic(modname, 'can''t get tnat',    1)
     422   IF(getKey('toce',    toce,    isoName)) CALL abort_physic(modname, 'can''t get toce',    1)
     423   IF(getKey('tcorr',   tcorr,   isoName)) CALL abort_physic(modname, 'can''t get tcorr',   1)
     424   IF(getKey('talph1',  talph1,  isoName)) CALL abort_physic(modname, 'can''t get talph1',  1)
     425   IF(getKey('talph2',  talph2,  isoName)) CALL abort_physic(modname, 'can''t get talph2',  1)
     426   IF(getKey('talph3',  talph3,  isoName)) CALL abort_physic(modname, 'can''t get talph3',  1)
     427   IF(getKey('talps1',  talps1,  isoName)) CALL abort_physic(modname, 'can''t get talps1',  1)
     428   IF(getKey('talps2',  talps2,  isoName)) CALL abort_physic(modname, 'can''t get talps2',  1)
     429   IF(getKey('tkcin0',  tkcin0,  isoName)) CALL abort_physic(modname, 'can''t get tkcin0',  1)
     430   IF(getKey('tkcin1',  tkcin1,  isoName)) CALL abort_physic(modname, 'can''t get tkcin1',  1)
     431   IF(getKey('tkcin2',  tkcin2,  isoName)) CALL abort_physic(modname, 'can''t get tkcin2',  1)
     432   IF(getKey('tdifrel', tdifrel, isoName)) CALL abort_physic(modname, 'can''t get tdifrel', 1)
     433   DO ixt = 1, niso
     434      IF     (ixt == iso_HTO) THEN; tdifrel(ixt) = 1./0.968
     435      ELSE IF(ixt == iso_HDO) THEN; tdifrel(ixt) = 1./0.9755
     436      ELSE IF(ixt == iso_O17) THEN; tdifrel(ixt) = 1./0.98555
     437      ELSE IF(ixt == iso_O18) THEN; tdifrel(ixt) = 1./0.9723
     438      ELSE IF(ixt == iso_eau) THEN; tdifrel(ixt) = 1.; END IF
     439   END DO
     440   IF(getKey('alpha_liq_sol', alpha_liq_sol, isoName)) CALL abort_physic(modname, 'can''t get alpha_liq_sol',  1)
     441   IF(getKey('Rdefault',Rdefault,isoName)) CALL abort_physic(modname, 'can''t get Rdefault',1)
     442   IF(getKey('Rmethox', Rmethox, isoName)) CALL abort_physic(modname, 'can''t get Rmethox', 1)
     443   IF(.NOT.Rdefault_smow) Rdefault(:) = 0.0
     444
     445   !===========================================================================================================================
     446   END IF
     447   !===========================================================================================================================
    415448
    416449   !--- Sensitivity test: no kinetic effect in sfc evaporation
    … …  
    423456   CALL msg('285: verif initialisation:', modname)
    424457   DO ixt=1,niso
    425       CALL msg(' * striso('//TRIM(sxt)//') = <'//TRIM(striso(ixt))//'>',   modname)
    426       CALL msg(  '   tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname)
    427 !     CALL msg('   alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname)
    428 !     CALL msg(       '   tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))),        modname)
    429 !     CALL msg(      '   tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))),       modname)
     458      sxt=int2str(ixt)
     459      CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>',  modname)
     460      CALL msg(  '    tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname)
     461!     CALL msg('    alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname)
     462!     CALL msg(        '   tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))),        modname)
     463!     CALL msg(       '   tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))),       modname)
    430464   END DO
    431465   CALL msg('69:     lambda = '//TRIM(real2str(lambda_sursat)), modname)
Note: See TracChangeset for help on using the changeset viewer.