Changeset 3852 for LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common
- Timestamp:
- Feb 22, 2021, 5:28:31 PM (3 years ago)
- Location:
- LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common
- Files:
-
- 6 edited
- 1 moved
-
infotrac.F90 (modified) (2 diffs)
-
initdynav.F90 (modified) (2 diffs)
-
inithist.F (modified) (2 diffs)
-
iso_verif_dyn.F90 (moved) (moved from LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/iso_verif_dyn.F) (1 diff)
-
massbarxy.F90 (modified) (1 diff)
-
writedynav.F90 (modified) (2 diffs)
-
writehist.F (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/infotrac.F90
r3851 r3852 1 ! $Id$2 !3 1 MODULE infotrac 4 2 5 ! nqtot : total number of tracers and higher order of moment, water vapor and liquid included 6 INTEGER, SAVE :: nqtot 7 !CR: on ajoute le nombre de traceurs de l eau 8 INTEGER, SAVE :: nqo 9 10 ! nbtr : number of tracers not including higher order of moment or water vapor or liquid 11 ! number of tracers used in the physics 12 INTEGER, SAVE :: nbtr 13 14 ! CRisi: nb traceurs pères= directement advectés par l'air 15 INTEGER, SAVE :: nqperes 16 17 ! Name variables 18 CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics 19 CHARACTER(len=23), ALLOCATABLE, DIMENSION(:), SAVE :: ttext ! tracer long name for diagnostics 20 21 ! iadv : index of trasport schema for each tracer 22 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iadv 23 24 ! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the 25 ! dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code. 26 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique 27 28 ! CRisi: tableaux de fils 29 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqfils 30 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les générations 31 INTEGER, SAVE :: nqdesc_tot 32 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqfils 33 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iqpere 34 REAL :: qperemin,masseqmin,ratiomin ! MVals et CRisi 35 PARAMETER (qperemin=1e-16,masseqmin=1e-16,ratiomin=1e-16) ! MVals 36 37 ! conv_flg(it)=0 : convection desactivated for tracer number it 38 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: conv_flg 39 ! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it 40 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: pbl_flg 41 42 CHARACTER(len=4),SAVE :: type_trac 43 CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym 44 45 ! CRisi: cas particulier des isotopes 46 LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso 47 INTEGER :: niso_possibles 48 PARAMETER ( niso_possibles=5) 49 REAL, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal 50 LOGICAL, DIMENSION(niso_possibles),SAVE :: use_iso 51 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase) 52 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot 53 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot 54 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numéro de la zone de tracage en fn de nqtot 55 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numéro de la zone de tracage en fn de nqtot 56 INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles 57 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: index_trac ! numéro ixt en fn izone, indnum entre 1 et niso 58 INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso 59 3 USE strings_mod, ONLY: msg, find, strIdx, strFind, strHead, dispTable, cat, get_in, & 4 fmsg, test, int2str, strParse, strTail, strReduce, strStack, modname, testFile 5 USE readTracFiles_mod, ONLY: readTracersFiles, getKey_init, nphases, delPhase, aliasTracer, & 6 tran0, readIsotopesFile, getKey, known_phases, addPhase, indexUpdate 7 USE trac_types_mod, ONLY: tra, iso, kys 8 9 IMPLICIT NONE 10 11 PRIVATE 12 13 !=== FOR TRACERS: 14 PUBLIC :: tra, tracers, type_trac !--- Derived type, full database, tracers type keyword 15 PUBLIC :: nqtot, nbtr, nqo !--- Main dimensions 16 PUBLIC :: infotrac_init, aliasTracer !--- Initialization, tracers alias creation 17 PUBLIC :: itr_indice !--- Indexes of the tracers passed to phytrac 18 PUBLIC :: niadv !--- Indexes of true tracers (<=nqtot, such that iadv(idx)>0) 19 PUBLIC :: solsym, conv_flg, pbl_flg 20 21 !=== FOR ISOTOPES: General 22 !--- General 23 PUBLIC :: iso, isotopes, nbIso !--- Derived type, full isotopes families database + nb of families 24 PUBLIC :: isoSelect , ixIso !--- Isotopes family selection tool + selected family index 25 !=== FOR ISOTOPES: Specific to H2O isotopes 26 PUBLIC :: iH2O, tnat, alpha_ideal !--- H2O isotopes index, natural abundance, fractionning coeff. 27 !=== FOR ISOTOPES: Depending on selected isotopes family 28 PUBLIC :: isotope, isoKeys !--- Selected isotopes database + associated keys (cf. getKey) 29 PUBLIC :: isoName, isoZone, isoPhas !--- Isotopes and tagging zones names, phases 30 PUBLIC :: niso, nzon, npha, nitr !--- " " numbers + isotopes & tagging tracers number 31 PUBLIC :: iZonIso, iTraPha !--- 2D index tables to get "iq" index 32 PUBLIC :: isoCheck !--- Run isotopes checking routines 33 34 !=== FOR BOTH TRACERS AND ISOTOPES 35 PUBLIC :: getKey !--- Get a key from "tracers" or "isotope" 36 37 !=== FOR STRATOSPHERIC AEROSOLS 60 38 #ifdef CPP_StratAer 61 !--CK/OB for stratospheric aerosols 62 INTEGER, SAVE :: nbtr_bin 63 INTEGER, SAVE :: nbtr_sulgas 64 INTEGER, SAVE :: id_OCS_strat 65 INTEGER, SAVE :: id_SO2_strat 66 INTEGER, SAVE :: id_H2SO4_strat 67 INTEGER, SAVE :: id_BIN01_strat 68 INTEGER, SAVE :: id_TEST_strat 69 #endif 70 39 PUBLIC :: nbtr_bin, nbtr_sulgas, id_OCS_strat, id_H2SO4_strat, id_SO2_strat, id_BIN01_strat, id_TEST_strat 40 #endif 41 42 INTERFACE isoSelect; MODULE PROCEDURE isoSelectByIndex, isoSelectByName; END INTERFACE isoSelect 43 44 !=== CONVENTIONS FOR TRACERS NUMBERS: 45 ! |--------------------+----------------------+-----------------+---------------+----------------------------| 46 ! | water in different | water tagging | water isotopes | other tracers | additional tracers moments | 47 ! | phases: H2O-[gls] | isotopes | | | for higher order schemes | 48 ! |--------------------+----------------------+-----------------+---------------+----------------------------| 49 ! | | | | | | 50 ! |<-- nqo -->|<-- nqo*niso* nzon -->|<-- nqo*niso -->|<-- nbtr -->|<-- (nmom) -->| 51 ! | | | | 52 ! | |<-- nqo*niso*(nzon+1) = nqo*nitr -->|<-- nqtottr = nbtr + nmom -->| 53 ! | = nqtot - nqo*(nitr+1) | 54 ! | | 55 ! |<-- nqtrue = nbtr + nqo*(nitr+1) -->| | 56 ! | | 57 ! |<-- nqtot = nqtrue + nmom -->| 58 ! | | 59 ! |----------------------------------------------------------------------------------------------------------| 60 ! NOTES FOR THIS TABLE: 61 ! * The used "niso", "nzon" and "nitr" are the H2O components of "isotopes(ip)" (isotopes(ip)%prnt == 'H2O'), 62 ! since water is so far the sole tracers family removed from the main tracers table. 63 ! * For water, "nqo" is equal to the more general field "isotopes(ip)%npha". 64 ! * "niso", "nzon", "nitr", "npha" are defined for other isotopic tracers families, if any. 65 ! 66 !=== TRACERS DESCRIPTOR: DERIVED TYPE EMBEDDING MOST OF THE USEFUL QUANTITIES (LENGTH: nqtot) 67 ! Each entry is accessible using "%" sign. 68 ! |------------+-------------------------------------------------+-------------+------------------------+ 69 ! | entry | Meaning | Former name | Possible values | 70 ! |------------+-------------------------------------------------+-------------+------------------------+ 71 ! | name | Name (short) | tname | | 72 ! | nam1 | Name of the 1st generation ancestor | / | | 73 ! | prnt | Name of the parent | / | | 74 ! | lnam | Long name (with adv. scheme suffix) for outputs | ttext | | 75 ! | type | Type (so far: tracer or tag) | / | tracer,tag | 76 ! | phas | Phases list ("g"as / "l"iquid / "s"olid) | / | [g][l][s] | 77 ! | iadv | Advection scheme number | iadv | 1-20,30 exc. 3-9,15,19 | 78 ! | igen | Generation (>=1) | / | | 79 ! | itr | Index in "tr_seri" (0: absent from physics) | cf. niadv | 1:nqtottr | 80 ! | iprnt | Index of the parent tracer | iqpere | 1:nqtot | 81 ! | idesc | Indexes of the childs (all generations) | iqfils | 1:nqtot | 82 ! | ndesc | Number of the descendants (all generations) | nqdesc | 1:nqtot | 83 ! | nchld | Number of childs (first generation only) | nqfils | 1:nqtot | 84 ! | keys | key/val pairs accessible with "getKey" routine | / | | 85 ! | iso_num | Isotope name index in iso(igr)%name(:) | iso_indnum | 1:niso | 86 ! | iso_zon | Isotope zone index in iso(igr)%zone(:) | zone_num | 1:nzon | 87 ! | iso_pha | Isotope phase index in iso(igr)%phas | phase_num | 1:npha | 88 ! +------------+-------------------------------------------------+-------------+------------------------+ 89 ! 90 !=== ISOTOPES DESCRIPTOR: DERIVED TYPE EMBEDDING MOST OF THE USEFUL QUANTITIES (LENGTH: NUMBER OF ISOTOPES FAMILIES USED) 91 ! Each entry is accessible using "%" sign. 92 ! |------------+-------------------------------------------------+-------------+-----------------------+ 93 ! | entry | Meaning | Former name | Possible values | 94 ! |------------+-------------------------------------------------+-------------+-----------------------+ 95 ! | prnt | Parent tracer (isotopes family name) | | | 96 ! | trac, nitr | Isotopes & tagging tracers + number of elements | | | 97 ! | zone, nzon | Geographic tagging zones + number of elements | | | 98 ! | phas, npha | Phases list + number of elements | | [g][l][s], 1:3 | 99 ! | niso | Number of isotopes, excluding tagging tracers | | | 100 ! | iTraPha | Index in "xt" = f(iname(niso+1:nitr),iphas) | iqiso | 1:niso | 101 ! | iZonIso | Index in "xt" = f(izone, iname(1:niso)) | index_trac | 1:nzon | 102 ! |------------+-------------------------------------------------+-------------+-----------------------+ 103 104 105 106 !=== DIMENSIONS OF THE TRACERS TABLES AND OTHER SCALAR VARIABLES 107 INTEGER, SAVE :: nqtot, & !--- Tracers nb in dynamics (incl. higher moments & water) 108 nbtr, & !--- Tracers nb in physics (excl. higher moments & water) 109 nqo, & !--- Number of water phases 110 nbIso !--- Number of available isotopes family 111 CHARACTER(LEN=256), SAVE :: type_trac !--- Keyword for tracers type 112 !$OMP THREADPRIVATE(nqtot, nbtr, nqo, nbIso, type_trac) 113 114 !=== DERIVED TYPES EMBEDDING MOST INFORMATIONS ABOUT TRACERS AND ISOTOPES FAMILIES 115 TYPE(tra), TARGET, SAVE, ALLOCATABLE :: tracers(:) !=== TRACERS DESCRIPTORS VECTOR 116 TYPE(iso), TARGET, SAVE, ALLOCATABLE :: isotopes(:) !=== ISOTOPES PARAMETERS VECTOR 117 !$OMP THREADPRIVATE(tracers, isotopes) 118 119 !=== ALIASES FOR CURRENTLY SELECTED ISOTOPES FAMILY OF VARIABLES EMBEDDED IN THE VECTOR "isotopes" 120 TYPE(iso), SAVE, POINTER :: isotope !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR 121 INTEGER, SAVE :: ixIso, iH2O !--- Index of the selected isotopes family and H2O family 122 LOGICAL, SAVE, POINTER :: isoCheck !--- Flag to trigger the checking routines 123 TYPE(kys), SAVE, POINTER :: isoKeys(:) !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName) 124 CHARACTER(LEN=256), SAVE, POINTER :: isoName(:), & !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY 125 isoZone(:), & !--- TAGGING ZONES FOR THE CURRENTLY SELECTED FAMILY 126 isoPhas !--- USED PHASES FOR THE CURRENTLY SELECTED FAMILY 127 INTEGER, SAVE :: niso, nzon, npha, & !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES 128 nitr !--- NUMBER OF ISOTOPES + ISOTOPIC TAGGING TRACERS 129 INTEGER, SAVE, POINTER :: iZonIso(:,:) !--- INDEX IN "isoTrac" AS f(tagging zone, isotope) 130 INTEGER, SAVE, POINTER :: iTraPha(:,:) !=== INDEX IN "isoTrac" AS f(isotopic tracer, phase) 131 !$OMP THREADPRIVATE(isotope, ixIso,iH2O, isoCheck, isoKeys, isoName,isoZone,isoPhas, niso,nzon,npha,nitr, iZonIso,iTraPha) 132 133 !=== VARIABLES EMBEDDED IN "tracers", BUT DUPLICATED, AS THEY ARE RATHER FREQUENTLY USED + VARIABLES FOR INCA 134 REAL, SAVE, ALLOCATABLE :: tnat(:), & !--- Natural relative abundance of water isotope (niso) 135 alpha_ideal(:) !--- Ideal fractionning coefficient (for initial state) (niso) 136 INTEGER, SAVE, ALLOCATABLE :: conv_flg(:), & !--- Convection activation ; needed for INCA (nbtr) 137 pbl_flg(:), & !--- Boundary layer activation ; needed for INCA (nbtr) 138 itr_indice(:), & !--- Indexes of the tracers passed to phytrac (nqtottr) 139 niadv(:) 140 CHARACTER(LEN=8), SAVE, ALLOCATABLE :: solsym(:) !--- Names from INCA (nbtr) 141 !OMP THREADPRIVATE(tnat, alpha_ideal, conv_flg, pbl_flg, itr_indice, solsym) 142 143 #ifdef CPP_StratAer 144 !=== SPECIFIC TO STRATOSPHERIC AEROSOLS (CK/OB) 145 INTEGER, SAVE :: nbtr_bin, nbtr_sulgas, id_OCS_strat, id_H2SO4_strat, id_SO2_strat, id_BIN01_strat, id_TEST_strat 146 !OMP THREADPRIVATE(nbtr_bin, nbtr_sulgas, id_OCS_strat, id_H2SO4_strat, id_SO2_strat, id_BIN01_strat, id_TEST_strat) 147 #endif 148 71 149 CONTAINS 72 150 73 74 151 SUBROUTINE infotrac_init 152 USE control_mod, ONLY: planet_type, config_inca 75 153 #ifdef REPROBUS 76 USE CHEM_REP, ONLY : Init_chem_rep_trac 77 #endif 78 IMPLICIT NONE 79 !======================================================================= 154 USE chem_rep, ONLY: Init_chem_rep_trac 155 #endif 156 !============================================================================================================================== 80 157 ! 81 158 ! Auteur: P. Le Van /L. Fairhead/F.Hourdin 82 159 ! ------- 83 ! Modif special traceur F.Forget 05/94 84 ! Modif M-A Filiberti 02/02 lecture de traceur.def 160 ! 161 ! Modifications: 162 ! -------------- 163 ! 05/94: F.Forget Modif special traceur 164 ! 02/02: M-A Filiberti Lecture de traceur.def 165 ! 06/20: D. Cugnet Nouveaux tracer.def et tracer_*.def + encapsulation (types tr et iso) 85 166 ! 86 167 ! Objet: … … 88 169 ! GCM LMD nouvelle grille 89 170 ! 90 !======================================================================= 171 !============================================================================================================================== 91 172 ! ... modification de l'integration de q ( 26/04/94 ) .... 92 !----------------------------------------------------------------------- 93 ! Declarations 94 95 INCLUDE "dimensions.h" 96 INCLUDE "iniprint.h" 97 173 !------------------------------------------------------------------------------------------------------------------------------ 174 ! Declarations: 175 ! INCLUDE "dimensions.h" 176 177 !------------------------------------------------------------------------------------------------------------------------------ 98 178 ! Local variables 99 INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv ! index of horizontal trasport schema 100 INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv ! index of vertical trasport schema 101 102 INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca ! index of horizontal trasport schema 103 INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca ! index of vertical trasport schema 104 105 CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_0 ! tracer short name 106 CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi 107 CHARACTER(len=3), DIMENSION(30) :: descrq 108 CHARACTER(len=1), DIMENSION(3) :: txts 109 CHARACTER(len=2), DIMENSION(9) :: txtp 110 CHARACTER(len=23) :: str1,str2 111 112 INTEGER :: nqtrue ! number of tracers read from tracer.def, without higer order of moment 113 INTEGER :: iq, new_iq, iiq, jq, ierr 114 INTEGER :: ifils,ipere,generation ! CRisi 115 LOGICAL :: continu,nouveau_traceurdef 116 INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def 117 CHARACTER(len=15) :: tchaine 118 119 character(len=*),parameter :: modname="infotrac_init" 120 !----------------------------------------------------------------------- 179 INTEGER, ALLOCATABLE :: hadv(:), hadv_inca(:), & !--- Horizontal/vertical transport scheme number 180 vadv(:), vadv_inca(:) !--- + specific INCA versions 181 CHARACTER(LEN=1) :: ph !--- Phase 182 CHARACTER(LEN=2) :: suff(9) !--- Suffixes for schemes of order 3 or 4 (Prather) 183 CHARACTER(LEN=3) :: descrq(30) !--- Advection scheme description 184 CHARACTER(LEN=4) :: oldH2O(3) !--- Old water names 185 CHARACTER(LEN=256) :: newH2O, iname, isoPhase !--- New water and isotope names, phases list 186 CHARACTER(LEN=256) :: msg1, msg2 !--- Strings for messages 187 CHARACTER(LEN=256), ALLOCATABLE, DIMENSION(:) :: & !--- Temporary storage 188 isoName, isoZone, tra0, zon0, tag0, n, p, z, str 189 INTEGER :: fType !--- Tracers description file type ; 0: none 190 !--- 1: "traceur.def" 2: "tracer.def" 3: "tracer_*.def" 191 INTEGER :: nqtrue !--- Tracers nb from tracer.def (no higher order moments) 192 INTEGER :: iad !--- Advection scheme 193 INTEGER :: iH2O !--- Index in "isotopes(:)" of H2O family 194 INTEGER :: ic,ip,iq,jq, it,nt, im,nm, ix, iz, niso, nzone, ntiso !--- Indexes and temporary variables 195 LOGICAL, ALLOCATABLE :: lisoGen2(:), & !--- Mask for second generation isotopes 196 lisoName(:), & !--- Mask for water isotopes 197 lisoZone(:), ll(:) !--- Mask for water isotopes tagging tracers 198 LOGICAL :: lerr 199 TYPE(tra), ALLOCATABLE, TARGET :: ttr(:) 200 TYPE(tra), POINTER :: t1, t(:) 201 TYPE(iso), POINTER :: s 202 !------------------------------------------------------------------------------------------------------------------------------ 121 203 ! Initialization : 122 ! 123 txts=(/'x','y','z'/) 124 txtp=(/'x ','y ','z ','xx','xy','xz','yy','yz','zz'/) 125 126 descrq(14)='VLH' 127 descrq(10)='VL1' 128 descrq(11)='VLP' 129 descrq(12)='FH1' 130 descrq(13)='FH2' 131 descrq(16)='PPM' 132 descrq(17)='PPS' 133 descrq(18)='PPP' 134 descrq(20)='SLP' 135 descrq(30)='PRA' 136 137 138 ! Coherence test between parameter type_trac, config_inca and preprocessing keys 139 IF (type_trac=='inca') THEN 140 WRITE(lunout,*) 'You have choosen to couple with INCA chemestry model : type_trac=', & 141 type_trac,' config_inca=',config_inca 142 IF (config_inca/='aero' .AND. config_inca/='aeNP' .AND. config_inca/='chem') THEN 143 WRITE(lunout,*) 'Incoherence between type_trac and config_inca. Model stops. Modify run.def' 144 CALL abort_gcm('infotrac_init','Incoherence between type_trac and config_inca',1) 145 END IF 204 !------------------------------------------------------------------------------------------------------------------------------ 205 modname = 'infotrac_init' 206 type_trac='lmdz'!'lmdz,inca' 207 suff = ['x ','y ','z ','xx','xy','xz','yy','yz','zz'] 208 descrq( 1: 2) = ['LMV','BAK'] 209 descrq(10:20) = ['VL1','VLP','FH1','FH2','VLH',' ','PPM','PPS','PPP',' ','SLP'] 210 descrq(30) = 'PRA' 211 oldH2O = ['H2Ov','H2Ol','H2Oi'] 212 213 !--- MESSAGE ABOUT THE CHOSEN CONFIGURATION 214 CALL msg('type_trac='//TRIM(type_trac)) 215 IF(strParse(type_trac, ',', str, n=nt)) CALL abort_gcm(modname,'can''t parse "type_trac = '//TRIM(type_trac)//'"',1) 216 DO it = 1, nt !--- nt>1 if "type_trac" is a coma-separated keywords list 217 msg1 = 'For type_trac = "'//TRIM(str(it))//'":' 218 SELECT CASE(str(it)) 219 CASE('inca'); CALL msg(TRIM(msg1)//' coupling with INCA chemistry model, config_inca='//config_inca) 220 CASE('repr'); CALL msg(TRIM(msg1)//' coupling with REPROBUS chemistry model') 221 CASE('co2i'); CALL msg(TRIM(msg1)//' you have chosen to run with CO2 cycle') 222 CASE('coag'); CALL msg(TRIM(msg1)//' tracers are treated for COAGULATION tests') 223 CASE('lmdz'); CALL msg(TRIM(msg1)//' tracers are treated in LMDZ only') 224 CASE DEFAULT 225 CALL abort_gcm(modname,'type_trac='//TRIM(str(it))//' not possible yet.',1) 226 END SELECT 227 END DO 228 229 !--- COHERENCE TEST BETWEEN "type_trac", "config_inca" AND PREPROCESSING KEYS 230 DO it=1,nt 231 SELECT CASE(type_trac) 232 CASE('inca'); IF(ALL(['aero', 'aeNP', 'chem']/=config_inca)) & 233 CALL abort_gcm(modname, 'Mismatch between type_trac and config_inca. Please modify "run.def"',1) 146 234 #ifndef INCA 147 WRITE(lunout,*) 'To run this option you must add cpp key INCA and compile with INCA code' 148 CALL abort_gcm('infotrac_init','You must compile with cpp key INCA',1) 149 #endif 150 ELSE IF (type_trac=='repr') THEN 151 WRITE(lunout,*) 'You have choosen to couple with REPROBUS chemestry model : type_trac=', type_trac 235 CALL abort_gcm(modname, 'You must add cpp key INCA and compile with INCA code',1) 236 #endif 237 CASE('repr') 152 238 #ifndef REPROBUS 153 WRITE(lunout,*) 'To run this option you must add cpp key REPROBUS and compile with REPRPBUS code' 154 CALL abort_gcm('infotrac_init','You must compile with cpp key REPROBUS',1) 155 #endif 156 ELSE IF (type_trac == 'co2i') THEN 157 WRITE(lunout,*) 'You have chosen to run with CO2 cycle: type_trac=', type_trac 158 ELSE IF (type_trac == 'coag') THEN 159 WRITE(lunout,*) 'Tracers are treated for COAGULATION tests : type_trac=', type_trac 239 CALL abort_gcm(modname, 'You must add cpp key REPROBUS and compile with REPROBUS code',1) 240 #endif 241 CASE('coag') 160 242 #ifndef CPP_StratAer 161 WRITE(lunout,*) 'To run this option you must add cpp key StratAer and compile with StratAer code' 162 CALL abort_gcm('infotrac_init','You must compile with cpp key StratAer',1) 163 #endif 164 ELSE IF (type_trac == 'lmdz') THEN 165 WRITE(lunout,*) 'Tracers are treated in LMDZ only : type_trac=', type_trac 166 ELSE 167 WRITE(lunout,*) 'type_trac=',type_trac,' not possible. Model stops' 168 CALL abort_gcm('infotrac_init','bad parameter',1) 169 END IF 170 171 ! Test if config_inca is other then none for run without INCA 172 IF (type_trac/='inca' .AND. config_inca/='none') THEN 173 WRITE(lunout,*) 'config_inca will now be changed to none as you do not couple with INCA model' 174 config_inca='none' 175 END IF 176 177 !----------------------------------------------------------------------- 178 ! 179 ! 1) Get the true number of tracers + water vapor/liquid 180 ! Here true tracers (nqtrue) means declared tracers (only first order) 181 ! 182 !----------------------------------------------------------------------- 183 IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN 184 OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr) 185 IF(ierr.EQ.0) THEN 186 WRITE(lunout,*) trim(modname),': Open traceur.def : ok' 187 READ(90,*) nqtrue 188 write(lunout,*) 'nqtrue=',nqtrue 189 ELSE 190 WRITE(lunout,*) trim(modname),': Problem in opening traceur.def' 191 WRITE(lunout,*) trim(modname),': WARNING using defaut values' 192 IF (planet_type=='earth') THEN 193 nqtrue=4 ! Default value for Earth 194 ELSE 195 nqtrue=1 ! Default value for other planets 196 ENDIF 197 ENDIF 198 !jyg< 199 !! if ( planet_type=='earth') then 200 !! ! For Earth, water vapour & liquid tracers are not in the physics 201 !! nbtr=nqtrue-2 202 !! else 203 !! ! Other planets (for now); we have the same number of tracers 204 !! ! in the dynamics than in the physics 205 !! nbtr=nqtrue 206 !! endif 207 !>jyg 208 ELSE ! type_trac=inca 209 !jyg< 210 ! The traceur.def file is used to define the number "nqo" of water phases 211 ! present in the simulation. Default : nqo = 2. 212 OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr) 213 IF(ierr.EQ.0) THEN 214 WRITE(lunout,*) trim(modname),': Open traceur.def : ok' 215 READ(90,*) nqo 216 ELSE 217 WRITE(lunout,*) trim(modname),': Using default value for nqo' 218 nqo=2 219 ENDIF 220 IF (nqo /= 2 .AND. nqo /= 3 ) THEN 221 WRITE(lunout,*) trim(modname),': nqo=',nqo, ' is not allowded. Only 2 or 3 water phases allowed' 222 CALL abort_gcm('infotrac_init','Bad number of water phases',1) 223 END IF 224 ! nbtr has been read from INCA by init_const_lmdz() in gcm.F 243 CALL abort_gcm(modname, 'You must add cpp key StratAer and compile with StratAer code',1) 244 #endif 245 END SELECT 246 END DO 247 248 !--- Disable "config_inca" option for a run without INCA if it differs from "none" 249 IF (ALL(str(:) /= 'inca') .AND. config_inca /= 'none') THEN 250 CALL msg('setting config_inca="none" as you do not couple with INCA model') 251 config_inca = 'none' 252 END IF 253 254 !------------------------------------------------------------------------------------------------------------------------------ 255 ! 1) Get the numbers of: true tracers "nqtrue", water tracers "nqo" (vapor/liquid/solid) 256 ! (here, "true" tracers means declared tracers, first order only) 257 ! Deal with the advection scheme choice for water and tracers: 258 ! iadv = 1 "LMDZ-specific humidity transport" (for H2O vapour) LMV 259 ! iadv = 2 backward (for H2O liquid) BAK 260 ! iadv = 14 Van-Leer + specific humidity, modified by Francis Codron VLH 261 ! iadv = 10 Van-Leer (chosen for vapour and liquid water) VL1 262 ! iadv = 11 Van-Leer for hadv and PPM version (Monotonic) for vadv VLP 263 ! iadv = 12 Frederic Hourdin I FH1 264 ! iadv = 13 Frederic Hourdin II FH2 265 ! iadv = 16 Monotonic PPM (Collela & Woodward 1984) PPM 266 ! iadv = 17 Semi-monotonic PPM (overshoots allowed) PPS 267 ! iadv = 18 Definite positive PPM (overshoots and undershoots allowed) PPP 268 ! iadv = 20 Slopes SLP 269 ! iadv = 30 Prather PRA 270 ! 271 ! In array q(ij,l,iq) : iq = 1 for vapour water 272 ! iq = 2 for liquid water 273 ! [iq = 3 for ice water] 274 ! And optionaly: iq = 3[4],nqtot for other tracers 275 !------------------------------------------------------------------------------------------------------------------------------ 276 ! Get choice of advection scheme from file tracer.def or from INCA 277 !------------------------------------------------------------------------------------------------------------------------------ 278 279 IF(readTracersFiles(type_trac, fType, tracers)) CALL abort_gcm(modname,'problem with tracers file(s)',1) 280 CALL msg(fType == 0, 'WARNING: USING DEFAULT VALUES !') 281 282 !---------------------------------------------------------------------------------------------------------------------------- 283 SELECT CASE(fType) 284 !---------------------------------------------------------------------------------------------------------------------------- 285 CASE(0) !=== NO READABLE TRACERS CONFIG FILE => DEFAULT 286 !-------------------------------------------------------------------------------------------------------------------------- 287 IF(planet_type=='earth') THEN !--- Default for Earth 288 nqo = 2; nbtr = 2 289 tracers(:)%name = ['H2O-g','H2O-l','RN ','PB '] 290 tracers(:)%prnt = [tran0 ,tran0 ,tran0 ,tran0 ] 291 tracers(:)%igen = [1 ,1 ,1 ,1 ] 292 hadv = [14 ,10 ,10 ,10 ] 293 vadv = [14 ,10 ,10 ,10 ] 294 ELSE !--- Default for other planets 295 nqo = 0; nbtr = 1 296 tracers(:)%name = ['dummy'] 297 tracers(:)%prnt = ['dummy'] 298 tracers(:)%igen = [1 ] 299 hadv = [10 ] 300 vadv = [10 ] 301 END IF 302 nqtrue = nbtr + nqo 303 !-------------------------------------------------------------------------------------------------------------------------- 304 CASE(1) 305 !-------------------------------------------------------------------------------------------------------------------------- 306 IF(type_trac=='inca') THEN !=== OLD STYLE "traceur.def" FOR INCA FOUND 307 !------------------------------------------------------------------------------------------------------------------------ 308 nqo = SIZE(tracers(:), DIM=1) 309 WRITE(msg1,'(a,i0)')'Only 2 or 3 water phases allowed ; found nqo=',nqo 310 IF(nqo/=2 .AND. nqo/=3) CALL abort_gcm(modname,TRIM(msg1),1) 225 311 #ifdef INCA 226 CALL Init_chem_inca_trac(nbtr) 227 #endif 228 nqtrue=nbtr+nqo 229 230 ALLOCATE(hadv_inca(nbtr), vadv_inca(nbtr)) 231 232 ENDIF ! type_trac 233 !>jyg 234 235 IF ((planet_type=="earth").and.(nqtrue < 2)) THEN 236 WRITE(lunout,*) trim(modname),': nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum' 237 CALL abort_gcm('infotrac_init','Not enough tracers',1) 238 END IF 239 240 !jyg< 241 ! Transfert number of tracers to Reprobus 242 !! IF (type_trac == 'repr') THEN 243 !!#ifdef REPROBUS 244 !! CALL Init_chem_rep_trac(nbtr) 245 !!#endif 246 !! END IF 247 !>jyg 248 249 ! 250 ! Allocate variables depending on nqtrue 251 ! 252 ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue)) 253 254 ! 255 !jyg< 256 !! ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr)) 257 !! conv_flg(:) = 1 ! convection activated for all tracers 258 !! pbl_flg(:) = 1 ! boundary layer activated for all tracers 259 !>jyg 260 261 !----------------------------------------------------------------------- 262 ! 2) Choix des schemas d'advection pour l'eau et les traceurs 263 ! 264 ! iadv = 1 schema transport type "humidite specifique LMD" 265 ! iadv = 2 schema amont 266 ! iadv = 14 schema Van-leer + humidite specifique 267 ! Modif F.Codron 268 ! iadv = 10 schema Van-leer (retenu pour l'eau vapeur et liquide) 269 ! iadv = 11 schema Van-Leer pour hadv et version PPM (Monotone) pour vadv 270 ! iadv = 12 schema Frederic Hourdin I 271 ! iadv = 13 schema Frederic Hourdin II 272 ! iadv = 16 schema PPM Monotone(Collela & Woodward 1984) 273 ! iadv = 17 schema PPM Semi Monotone (overshoots autorisés) 274 ! iadv = 18 schema PPM Positif Defini (overshoots undershoots autorisés) 275 ! iadv = 20 schema Slopes 276 ! iadv = 30 schema Prather 277 ! 278 ! Dans le tableau q(ij,l,iq) : iq = 1 pour l'eau vapeur 279 ! iq = 2 pour l'eau liquide 280 ! Et eventuellement iq = 3,nqtot pour les autres traceurs 281 ! 282 ! iadv(1): choix pour l'eau vap. et iadv(2) : choix pour l'eau liq. 283 !------------------------------------------------------------------------ 284 ! 285 ! Get choice of advection schema from file tracer.def or from INCA 286 !--------------------------------------------------------------------- 287 IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN 288 IF(ierr.EQ.0) THEN 289 ! Continue to read tracer.def 290 DO iq=1,nqtrue 291 292 write(*,*) 'infotrac 237: iq=',iq 293 ! CRisi: ajout du nom du fluide transporteur 294 ! mais rester retro compatible 295 READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine 296 write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq) 297 write(lunout,*) 'tchaine=',trim(tchaine) 298 write(*,*) 'infotrac 238: IOstatus=',IOstatus 299 if (IOstatus.ne.0) then 300 CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1) 301 endif 302 ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un 303 ! espace ou pas au milieu de la chaine. 304 continu=.true. 305 nouveau_traceurdef=.false. 306 iiq=1 307 do while (continu) 308 if (tchaine(iiq:iiq).eq.' ') then 309 nouveau_traceurdef=.true. 310 continu=.false. 311 else if (iiq.lt.LEN_TRIM(tchaine)) then 312 iiq=iiq+1 313 else 314 continu=.false. 315 endif 316 enddo 317 write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef 318 if (nouveau_traceurdef) then 319 write(lunout,*) 'C''est la nouvelle version de traceur.def' 320 tnom_0(iq)=tchaine(1:iiq-1) 321 tnom_transp(iq)=tchaine(iiq+1:15) 322 else 323 write(lunout,*) 'C''est l''ancienne version de traceur.def' 324 write(lunout,*) 'On suppose que les traceurs sont tous d''air' 325 tnom_0(iq)=tchaine 326 tnom_transp(iq) = 'air' 327 endif 328 write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>' 329 write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>' 330 331 END DO !DO iq=1,nqtrue 332 CLOSE(90) 333 334 ELSE ! Without tracer.def, set default values 335 if (planet_type=="earth") then 336 ! for Earth, default is to have 4 tracers 337 hadv(1) = 14 338 vadv(1) = 14 339 tnom_0(1) = 'H2Ov' 340 tnom_transp(1) = 'air' 341 hadv(2) = 10 342 vadv(2) = 10 343 tnom_0(2) = 'H2Ol' 344 tnom_transp(2) = 'air' 345 hadv(3) = 10 346 vadv(3) = 10 347 tnom_0(3) = 'RN' 348 tnom_transp(3) = 'air' 349 hadv(4) = 10 350 vadv(4) = 10 351 tnom_0(4) = 'PB' 352 tnom_transp(4) = 'air' 353 else ! default for other planets 354 hadv(1) = 10 355 vadv(1) = 10 356 tnom_0(1) = 'dummy' 357 tnom_transp(1) = 'dummy' 358 endif ! of if (planet_type=="earth") 359 END IF 360 361 WRITE(lunout,*) trim(modname),': Valeur de traceur.def :' 362 WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue 363 DO iq=1,nqtrue 364 WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq) 365 END DO 366 367 IF ( planet_type=='earth') THEN 368 !CR: nombre de traceurs de l eau 369 IF (tnom_0(3) == 'H2Oi') THEN 370 nqo=3 371 ELSE 372 nqo=2 373 ENDIF 374 ! For Earth, water vapour & liquid tracers are not in the physics 375 nbtr=nqtrue-nqo 376 ELSE 377 ! Other planets (for now); we have the same number of tracers 378 ! in the dynamics than in the physics 379 nbtr=nqtrue 380 ENDIF 312 CALL Init_chem_inca_trac(nbtr) !--- Get nbtr from INCA 313 #endif 314 ALLOCATE(hadv_inca(nbtr), vadv_inca(nbtr), conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr)) 315 #ifdef INCA 316 !--- Activation of: Convection, Boundary layer 317 CALL init_transport(hadv_inca, vadv_inca, conv_flg, pbl_flg, solsym) 318 #endif 319 nqtrue = nbtr + nqo !--- Total number of tracers 320 ALLOCATE(ttr(nqtrue)); ttr(1:nqo) = tracers(1:nqo) 321 DO iq = nqo+1, nqtrue 322 ttr(iq)%name = solsym(iq) 323 ttr(iq)%prnt = tran0 324 ttr(iq)%igen = 1 325 hadv = hadv_inca(iq-nqo) 326 vadv = vadv_inca(iq-nqo) 327 END DO 328 CALL MOVE_ALLOC(FROM=ttr, TO=tracers) 329 !------------------------------------------------------------------------------------------------------------------------ 330 ELSE !=== OLD STYLE "traceur.def" CONFIG FILE FOUND 331 !------------------------------------------------------------------------------------------------------------------------ 332 nqo = 0 333 DO ip = 1, SIZE(oldH2O) 334 ix = strIdx(tracers(:)%name,oldH2O(ip)) !--- Old name of water in a specific phase (ix/=0) 335 IF(ix == 0) CYCLE 336 newH2O = 'H2O-'//known_phases(ip:ip) !--- Corresponding new name 337 nqo = nqo+1; tracers(ix)%name = newH2O !--- One more water phase ; replace old name with one 338 tracers(strFind(tracers(:)%nam1,oldH2O(ip)))%nam1 = newH2O 339 tracers(strFind(tracers(:)%prnt,oldH2O(ip)))%prnt = newH2O 340 END DO 341 nqtrue = SIZE(tracers,DIM=1) 342 nbtr = nqtrue - nqo 343 END IF 344 !-------------------------------------------------------------------------------------------------------------------------- 345 CASE DEFAULT !=== FOUND NEW STYLE TRACERS CONFIG FILE(S) 346 !-------------------------------------------------------------------------------------------------------------------------- 347 nqo = 2; IF(ANY(tracers(:)%name == 'H2O-s')) nqo=3 348 nqtrue = SIZE(tracers, DIM=1) 349 nbtr = nqtrue - nqo 350 !---------------------------------------------------------------------------------------------------------------------------- 351 END SELECT 352 !---------------------------------------------------------------------------------------------------------------------------- 353 CALL getKey_init(tracers) 354 IF(.NOT.ALLOCATED(hadv)) lerr = getKey('hadv', hadv) 355 IF(.NOT.ALLOCATED(vadv)) lerr = getKey('vadv', vadv) 356 IF(.NOT.ALLOCATED(solsym)) ALLOCATE(solsym(nbtr)) 357 IF(.NOT.ALLOCATED(conv_flg)) conv_flg = [(1, it=1, nbtr)] 358 IF(.NOT.ALLOCATED( pbl_flg)) pbl_flg = [(1, it=1, nbtr)] 381 359 382 360 #ifdef CPP_StratAer 383 IF (type_trac == 'coag') THEN 384 nbtr_bin=0 385 nbtr_sulgas=0 386 DO iq=1,nqtrue 387 IF (tnom_0(iq)(1:3)=='BIN') THEN !check if tracer name contains 'BIN' 388 nbtr_bin=nbtr_bin+1 389 ENDIF 390 IF (tnom_0(iq)(1:3)=='GAS') THEN !check if tracer name contains 'GAS' 391 nbtr_sulgas=nbtr_sulgas+1 392 ENDIF 393 ENDDO 394 print*,'nbtr_bin=',nbtr_bin 395 print*,'nbtr_sulgas=',nbtr_sulgas 396 DO iq=1,nqtrue 397 IF (tnom_0(iq)=='GASOCS') THEN 398 id_OCS_strat=iq-nqo 399 ENDIF 400 IF (tnom_0(iq)=='GASSO2') THEN 401 id_SO2_strat=iq-nqo 402 ENDIF 403 IF (tnom_0(iq)=='GASH2SO4') THEN 404 id_H2SO4_strat=iq-nqo 405 ENDIF 406 IF (tnom_0(iq)=='BIN01') THEN 407 id_BIN01_strat=iq-nqo 408 ENDIF 409 IF (tnom_0(iq)=='GASTEST') THEN 410 id_TEST_strat=iq-nqo 411 ENDIF 412 ENDDO 413 print*,'id_OCS_strat =',id_OCS_strat 414 print*,'id_SO2_strat =',id_SO2_strat 415 print*,'id_H2SO4_strat=',id_H2SO4_strat 416 print*,'id_BIN01_strat=',id_BIN01_strat 417 ENDIF 418 #endif 419 420 ENDIF ! (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac = 'coag') 421 !jyg< 422 ! 423 ! Transfert number of tracers to Reprobus 424 IF (type_trac == 'repr') THEN 361 IF (type_trac == 'coag') THEN 362 nbtr_bin=0 363 nbtr_sulgas=0 364 DO iq = 1, nqtrue 365 IF(tracers(iq)%name(1:3)=='BIN') nbtr_bin = nbtr_bin +1 366 IF(tracers(iq)%name(1:3)=='GAS') nbtr_sulgas = nbtr_sulgas+1 367 SELECT CASE(tracers(iq)%name) 368 CASE('BIN01'); id_BIN01_strat = iq - nqo; CALL msg('id_BIN01_strat=', id_BIN01_strat) 369 CASE('GASOCS'); id_OCS_strat = iq - nqo; CALL msg('id_OCS_strat =', id_OCS_strat) 370 CASE('GASSO2'); id_SO2_strat = iq - nqo; CALL msg('id_SO2_strat =', id_SO2_strat) 371 CASE('GASH2SO4'); id_H2SO4_strat = iq - nqo; CALL msg('id_H2SO4_strat=', id_H2SO4_strat) 372 CASE('GASTEST'); id_TEST_strat = iq - nqo; CALL msg('id_TEST_strat=' , id_TEST_strat) 373 END SELECT 374 END DO 375 CALL msg('nbtr_bin =',nbtr_bin) 376 CALL msg('nbtr_sulgas =',nbtr_sulgas) 377 END IF 378 #endif 379 380 !--- Transfert number of tracers to Reprobus 425 381 #ifdef REPROBUS 426 CALL Init_chem_rep_trac(nbtr,nqo,tnom_0) 427 #endif 428 END IF 429 ! 430 ! Allocate variables depending on nbtr 431 ! 432 ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr)) 433 conv_flg(:) = 1 ! convection activated for all tracers 434 pbl_flg(:) = 1 ! boundary layer activated for all tracers 435 ! 436 !! ELSE ! type_trac=inca : config_inca='aero' ou 'chem' 437 ! 438 IF (type_trac == 'inca') THEN ! config_inca='aero' ou 'chem' 439 !>jyg 440 ! le module de chimie fournit les noms des traceurs 441 ! et les schemas d'advection associes. excepte pour ceux lus 442 ! dans traceur.def 443 IF (ierr .eq. 0) then 444 DO iq=1,nqo 445 446 write(*,*) 'infotrac 237: iq=',iq 447 ! CRisi: ajout du nom du fluide transporteur 448 ! mais rester retro compatible 449 READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine 450 write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq) 451 write(lunout,*) 'tchaine=',trim(tchaine) 452 write(*,*) 'infotrac 238: IOstatus=',IOstatus 453 if (IOstatus.ne.0) then 454 CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1) 455 endif 456 ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un 457 ! espace ou pas au milieu de la chaine. 458 continu=.true. 459 nouveau_traceurdef=.false. 460 iiq=1 461 do while (continu) 462 if (tchaine(iiq:iiq).eq.' ') then 463 nouveau_traceurdef=.true. 464 continu=.false. 465 else if (iiq.lt.LEN_TRIM(tchaine)) then 466 iiq=iiq+1 467 else 468 continu=.false. 469 endif 470 enddo 471 write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef 472 if (nouveau_traceurdef) then 473 write(lunout,*) 'C''est la nouvelle version de traceur.def' 474 tnom_0(iq)=tchaine(1:iiq-1) 475 tnom_transp(iq)=tchaine(iiq+1:15) 476 else 477 write(lunout,*) 'C''est l''ancienne version de traceur.def' 478 write(lunout,*) 'On suppose que les traceurs sont tous d''air' 479 tnom_0(iq)=tchaine 480 tnom_transp(iq) = 'air' 481 endif 482 write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>' 483 write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>' 484 485 END DO !DO iq=1,nqtrue 486 CLOSE(90) 487 ELSE !! if traceur.def doesn't exist 488 tnom_0(1)='H2Ov' 489 tnom_transp(1) = 'air' 490 tnom_0(2)='H2Ol' 491 tnom_transp(2) = 'air' 492 hadv(1) = 10 493 hadv(2) = 10 494 vadv(1) = 10 495 vadv(2) = 10 496 ENDIF 497 498 #ifdef INCA 499 CALL init_transport( & 500 hadv_inca, & 501 vadv_inca, & 502 conv_flg, & 503 pbl_flg, & 504 solsym) 505 #endif 506 507 508 !jyg< 509 DO iq = nqo+1, nqtrue 510 hadv(iq) = hadv_inca(iq-nqo) 511 vadv(iq) = vadv_inca(iq-nqo) 512 tnom_0(iq)=solsym(iq-nqo) 513 tnom_transp(iq) = 'air' 514 END DO 515 516 END IF ! (type_trac == 'inca') 517 518 !----------------------------------------------------------------------- 519 ! 520 ! 3) Verify if advection schema 20 or 30 choosen 382 IF(type_trac == 'repr') CALL Init_chem_rep_trac(nbtr,nqo,tracers(:)%name) 383 #endif 384 385 !------------------------------------------------------------------------------------------------------------------------------ 386 ! 2) Verify if the advection scheme 20 or 30 have been chosen. 521 387 ! Calculate total number of tracers needed: nqtot 522 388 ! Allocate variables depending on total number of tracers 523 !----------------------------------------------------------------------- 524 new_iq=0 525 DO iq=1,nqtrue 526 ! Add tracers for certain advection schema 527 IF (hadv(iq)<20 .AND. vadv(iq)<20 ) THEN 528 new_iq=new_iq+1 ! no tracers added 529 ELSE IF (hadv(iq)==20 .AND. vadv(iq)==20 ) THEN 530 new_iq=new_iq+4 ! 3 tracers added 531 ELSE IF (hadv(iq)==30 .AND. vadv(iq)==30 ) THEN 532 new_iq=new_iq+10 ! 9 tracers added 533 ELSE 534 WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq) 535 CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1) 536 END IF 389 !------------------------------------------------------------------------------------------------------------------------------ 390 DO iq = 1, nqtrue 391 t1 => tracers(iq) 392 IF( hadv(iq)<20 .OR. (ANY(hadv(iq)==[20,30]) .AND. hadv(iq)==vadv(iq)) ) CYCLE 393 WRITE(msg1,'(2(a,i0))')' is not available: hadv=',hadv(iq),', vadv=',vadv(iq) 394 CALL msg('This choice of advection scheme for "'//TRIM(t1%name)//'"'//TRIM(msg1)) 395 CALL abort_gcm(modname,'Bad choice of advection scheme',1) 396 END DO 397 nqtot = COUNT( hadv< 20 .AND. vadv< 20 ) & !--- No additional tracer 398 + 4*COUNT( hadv==20 .AND. vadv==20 ) & !--- 3 additional tracers 399 + 10*COUNT( hadv==30 .AND. vadv==30 ) !--- 9 additional tracers 400 401 ! More tracers due to the choice of advection scheme => assign total number of tracers 402 IF( nqtot /= nqtrue ) THEN 403 CALL msg('The choice of advection scheme for one or more tracers makes it necessary to add tracers') 404 CALL msg('The number of true tracers is '//TRIM(int2str(nqtrue))) 405 CALL msg('The total number of tracers needed is '//TRIM(int2str(nqtot))) 406 END IF 407 ALLOCATE(ttr(nqtot)) 408 409 !------------------------------------------------------------------------------------------------------------------------------ 410 ! 3) Determine iadv, long and short name, generation number, phase and region 411 !------------------------------------------------------------------------------------------------------------------------------ 412 jq = 0; ttr(:)%iadv = -1 413 DO iq = 1, nqtrue 414 jq = jq + 1 415 t1 => tracers(iq) 416 417 !--- Verify choice of advection schema 418 iad = -1 419 IF(hadv(iq) == vadv(iq) ) iad = hadv(iq) 420 IF(hadv(iq)==10 .AND. vadv(iq)==16) iad = 11 421 CALL msg(iad == -1, 'This choice of advection scheme for "'//TRIM(t1%name)//'" '//'is not available: hadv = ' & 422 //TRIM(int2str(hadv(iq)))//', vadv='//TRIM(int2str(vadv(iq))) ) 423 IF(iad == -1) CALL abort_gcm(modname,'Bad choice of advection scheme - 2',1) 424 t1%lnam = t1%name; IF(iad /= 0) t1%lnam=TRIM(t1%name)//descrq(iad) 425 426 !--- Defining most fields of the tracer derived type 427 ttr(jq)%name = t1%name 428 ttr(jq)%nam1 = t1%nam1 429 ttr(jq)%prnt = t1%prnt 430 ttr(jq)%lnam = t1%lnam 431 ttr(jq)%type = t1%type 432 ttr(jq)%phas = t1%phas 433 ttr(jq)%iadv = iad 434 ttr(jq)%igen = t1%igen 435 436 IF(ALL([20,30] /= iad)) CYCLE !--- 1st order scheme: finished 437 IF(iad == 20) nm = 3 !--- 2nd order scheme 438 IF(iad == 30) nm = 9 !--- 3rd order scheme 439 ttr(jq+1:jq+nm)%name = [ (TRIM(t1%name)//'-'//TRIM(suff(im)), im=1, nm) ] 440 ttr(jq+1:jq+nm)%nam1 = [ (TRIM(t1%nam1)//'-'//TRIM(suff(im)), im=1, nm) ] 441 ttr(jq+1:jq+nm)%lnam = [ (TRIM(t1%lnam)//'-'//TRIM(suff(im)), im=1, nm) ] 442 ttr(jq+1:jq+nm)%prnt = t1%prnt 443 ttr(jq+1:jq+nm)%type = t1%type 444 ttr(jq+1:jq+nm)%phas = t1%phas 445 ttr(jq+1:jq+nm)%iadv = -iad 446 ttr(jq+1:jq+nm)%igen = t1%igen 447 jq = jq + nm 448 END DO 449 DEALLOCATE(hadv, vadv) 450 451 !--- Determine parent and childs indexes 452 CALL indexUpdate(ttr) 453 454 !=== TEST ADVECTION SCHEME 455 DO iq=1,nqtot ; t1 => ttr(iq); iad = t1%iadv 456 WRITE(msg1,'(a,i0)')'This LMDZ version has not been tested for option iadv=',iad 457 WRITE(msg2,'(a,i2,a)')'iadv=',iad,' not implemented yet for' 458 459 !--- ONLY TESTED VALUES FOR TRACERS FOR NOW: iadv = 14, 10 (and 0) 460 IF(ALL( [10,14,0] /= iad) ) CALL abort_gcm(modname, TRIM(msg1)//' ; only iadv=10 and iadv=14 are tested !', 1) 461 462 !--- ONLY TESTED VALUES FOR CHILDS FOR NOW: iadv = 10 (CHILDS: TRACERS OF GENERATION GREATER THAN 1) 463 IF(fmsg(iad/=10.AND.t1%igen>1,'WARNING ! '//TRIM(msg2)//' childs. Setting iadv=10 for "'//TRIM(t1%name)//'".')) t1%iadv=10 464 465 !--- ONLY TESTED VALUES FOR PARENTS FOR NOW: iadv = 14, 10 (PARENTS: TRACERS OF GENERATION 1) 466 IF(t1%igen==1 .AND. ALL([10,14]/=iad)) CALL abort_gcm(modname, TRIM(msg2)//' parents: schemes 10 or 14 only !', 1) 467 468 !--- iadv = 14 IS ONLY VALID FOR WATER VAPOUR 469 IF(fmsg(iad==14 .AND. t1%name(1:5)/='H2O-g', 'WARNING ! '//TRIM(msg1)//', found for "' & 470 //TRIM(t1%name)//'" but only valid for water vapour ! Setting iadv=10 for "'//TRIM(t1%name)//'".')) t1%iadv=10 471 END DO 472 473 !=== DISPLAY THE RESULTING LIST 474 CALL msg('Information stored in infotrac :') 475 IF(dispTable('isssiii', ['iq ','name ','long name','parent ','iadv ','ipar ','igen '], & 476 cat(ttr(:)%name, ttr(:)%lnam, ttr(:)%prnt), cat([(iq, iq=1, nqtot)], ttr(:)%iadv, ttr(:)%iprnt, ttr(:)%igen))) & 477 CALL abort_gcm(modname,"problem with the tracers table content",1) 478 479 CALL MOVE_ALLOC(FROM=ttr, TO=tracers) 480 t => tracers 481 482 !=== VARIABLES RELATED TO GENERATIONS 483 niadv = PACK( [(iq,iq=1,nqtot)], MASK=t(:)%iadv>=0) !--- Indexes of "true" tracers 484 485 p = PACK(delPhase(t%prnt),MASK=t%type=='tracer'.AND.t%igen==2)!--- Parents of 2nd generation isotopes 486 CALL strReduce(p, nbIso) 487 ALLOCATE(isotopes(nbIso)) 488 489 IF(nbIso==0) RETURN !=== NO ISOTOPES: FINISHED 490 491 CALL msg('Isotopes families required: '//strStack(p)) 492 493 !--- ISOTOPES RELATED VARIABLES ; NULL OR EMPTY IF NO ISOTOPES 494 isotopes(:)%prnt = p 495 DO ip = 1, SIZE(p) !--- Loop on isotopes categories 496 s => isotopes(ip) 497 iname = s%prnt 498 499 !=== Geographic tagging tracers descending on tracer "iname": mask, names, number 500 lisoZone = t(:)%type=='tag' .AND. delPhase(t(:)%nam1) == iname .AND. t(:)%igen == 3 501 s%zone = PACK(strTail(t(:)%name,'_'), MASK = lisoZone) !--- Tagging zones names for isotopes category "iname" 502 CALL strReduce(s%zone) 503 s%nzon = SIZE(s%zone) !--- Tagging zones number for isotopes category "iname" 504 505 !=== Isotopes childs of tracer "iname": mask, names, number (same for each phase of "iname") 506 lisoName = t(:)%type=='tracer' .AND. delPhase(t(:)%prnt) == iname .AND. t(:)%phas == 'g' 507 ALLOCATE(s%keys(COUNT(lisoName))) 508 s%keys(:)%name = PACK(delPhase(t(:)%name), MASK = lisoName) !--- Effectively found isotopes of "iname" 509 s%niso = SIZE(s%keys) !--- Number of "effectively found isotopes of "iname" 510 s%trac = [s%keys%name, ((TRIM(s%keys(it)%name)//'_'//TRIM(s%zone(iz)), it=1, s%niso), iz=1, s%nzon)] 511 s%nitr = SIZE(s%trac) !--- " + their geographic tracers [ntraciso] 512 513 !=== Phases for tracer "iname" 514 s%phas = '' 515 DO ix = 1, nphases; IF(strIdx(t%name,addPhase(iname, known_phases(ix:ix))) /= 0) s%phas = TRIM(s%phas)//ph; END DO 516 s%npha = LEN_TRIM(s%phas) !--- Equal to "nqo" for water 517 518 !=== Tables giving the index in a table of effectively found items for each dynamical tracer (1<=iq<=nqtot) 519 DO iq = 1, nqtot 520 t1 => tracers(iq) 521 IF(t1%nam1 /= iname) CYCLE !--- Only deal with tracers descending on "iname" 522 t1%iso_igr = ip !--- Index of isotopes family in list "isotopes(:)%prnt" 523 t1%iso_num = strIdx(s%trac, delPhase(strHead(t1%name,'_')))!--- Index of current isotope in effective isotopes list 524 t1%iso_zon = strIdx(s%zone, strTail(t1%name,'_') )!--- Index of current isotope zone in effective zones list 525 t1%iso_pha = INDEX(s%phas,TRIM(t1%phas)) !--- Index of current isotope phase in effective phases list 526 IF(t1%igen /= 3) t1%iso_zon = 0 !--- Skip possible generation 2 tagging tracers 537 527 END DO 538 539 IF (new_iq /= nqtrue) THEN 540 ! The choice of advection schema imposes more tracers 541 ! Assigne total number of tracers 542 nqtot = new_iq 543 544 WRITE(lunout,*) trim(modname),': The choice of advection schema for one or more tracers' 545 WRITE(lunout,*) 'makes it necessary to add tracers' 546 WRITE(lunout,*) trim(modname)//': ',nqtrue,' is the number of true tracers' 547 WRITE(lunout,*) trim(modname)//': ',nqtot, ' is the total number of tracers needed' 548 549 ELSE 550 ! The true number of tracers is also the total number 551 nqtot = nqtrue 552 END IF 553 554 ! 555 ! Allocate variables with total number of tracers, nqtot 556 ! 557 ALLOCATE(tname(nqtot), ttext(nqtot)) 558 ALLOCATE(iadv(nqtot), niadv(nqtot)) 559 560 !----------------------------------------------------------------------- 561 ! 562 ! 4) Determine iadv, long and short name 563 ! 564 !----------------------------------------------------------------------- 565 new_iq=0 566 DO iq=1,nqtrue 567 new_iq=new_iq+1 568 569 ! Verify choice of advection schema 570 IF (hadv(iq)==vadv(iq)) THEN 571 iadv(new_iq)=hadv(iq) 572 ELSE IF (hadv(iq)==10 .AND. vadv(iq)==16) THEN 573 iadv(new_iq)=11 574 ELSE 575 WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq) 576 577 CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1) 578 END IF 579 580 str1=tnom_0(iq) 581 tname(new_iq)= tnom_0(iq) 582 IF (iadv(new_iq)==0) THEN 583 ttext(new_iq)=trim(str1) 584 ELSE 585 ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq)) 586 END IF 587 588 ! schemas tenant compte des moments d'ordre superieur 589 str2=ttext(new_iq) 590 IF (iadv(new_iq)==20) THEN 591 DO jq=1,3 592 new_iq=new_iq+1 593 iadv(new_iq)=-20 594 ttext(new_iq)=trim(str2)//txts(jq) 595 tname(new_iq)=trim(str1)//txts(jq) 596 END DO 597 ELSE IF (iadv(new_iq)==30) THEN 598 DO jq=1,9 599 new_iq=new_iq+1 600 iadv(new_iq)=-30 601 ttext(new_iq)=trim(str2)//txtp(jq) 602 tname(new_iq)=trim(str1)//txtp(jq) 603 END DO 604 END IF 605 END DO 606 607 ! 608 ! Find vector keeping the correspodence between true and total tracers 609 ! 610 niadv(:)=0 611 iiq=0 612 DO iq=1,nqtot 613 IF(iadv(iq).GE.0) THEN 614 ! True tracer 615 iiq=iiq+1 616 niadv(iiq)=iq 617 ENDIF 618 END DO 619 620 621 WRITE(lunout,*) trim(modname),': Information stored in infotrac :' 622 WRITE(lunout,*) trim(modname),': iadv niadv tname ttext :' 623 DO iq=1,nqtot 624 WRITE(lunout,*) iadv(iq),niadv(iq),& 625 ' ',trim(tname(iq)),' ',trim(ttext(iq)) 626 END DO 627 628 ! 629 ! Test for advection schema. 630 ! This version of LMDZ only garantees iadv=10 and iadv=14 (14 only for water vapour) . 631 ! 632 DO iq=1,nqtot 633 IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN 634 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' 635 CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1) 636 ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN 637 WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' 638 CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1) 639 END IF 640 END DO 641 642 643 ! CRisi: quels sont les traceurs fils et les traceurs pères. 644 ! initialiser tous les tableaux d'indices liés aux traceurs familiaux 645 ! + vérifier que tous les pères sont écrits en premières positions 646 ALLOCATE(nqfils(nqtot),nqdesc(nqtot)) 647 ALLOCATE(iqfils(nqtot,nqtot)) 648 ALLOCATE(iqpere(nqtot)) 649 nqperes=0 650 nqfils(:)=0 651 nqdesc(:)=0 652 iqfils(:,:)=0 653 iqpere(:)=0 654 nqdesc_tot=0 655 DO iq=1,nqtot 656 if (tnom_transp(iq) == 'air') then 657 ! ceci est un traceur père 658 WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere' 659 nqperes=nqperes+1 660 iqpere(iq)=0 661 else !if (tnom_transp(iq) == 'air') then 662 ! ceci est un fils. Qui est son père? 663 WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un fils' 664 continu=.true. 665 ipere=1 666 do while (continu) 667 if (tnom_transp(iq) == tnom_0(ipere)) then 668 ! Son père est ipere 669 WRITE(lunout,*) 'Le traceur',iq,'appele ', & 670 & trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere)) 671 nqfils(ipere)=nqfils(ipere)+1 672 iqfils(nqfils(ipere),ipere)=iq 673 iqpere(iq)=ipere 674 continu=.false. 675 else !if (tnom_transp(iq) == tnom_0(ipere)) then 676 ipere=ipere+1 677 if (ipere.gt.nqtot) then 678 WRITE(lunout,*) 'Le traceur',iq,'appele ', & 679 & trim(tnom_0(iq)),', est orphelin.' 680 CALL abort_gcm('infotrac_init','Un traceur est orphelin',1) 681 endif !if (ipere.gt.nqtot) then 682 endif !if (tnom_transp(iq) == tnom_0(ipere)) then 683 enddo !do while (continu) 684 endif !if (tnom_transp(iq) == 'air') then 685 enddo !DO iq=1,nqtot 686 WRITE(lunout,*) 'infotrac: nqperes=',nqperes 687 WRITE(lunout,*) 'nqfils=',nqfils 688 WRITE(lunout,*) 'iqpere=',iqpere 689 WRITE(lunout,*) 'iqfils=',iqfils 690 691 ! Calculer le nombre de descendants à partir de iqfils et de nbfils 692 DO iq=1,nqtot 693 generation=0 694 continu=.true. 695 ifils=iq 696 do while (continu) 697 ipere=iqpere(ifils) 698 if (ipere.gt.0) then 699 nqdesc(ipere)=nqdesc(ipere)+1 700 nqdesc_tot=nqdesc_tot+1 701 iqfils(nqdesc(ipere),ipere)=iq 702 ifils=ipere 703 generation=generation+1 704 else !if (ipere.gt.0) then 705 continu=.false. 706 endif !if (ipere.gt.0) then 707 enddo !do while (continu) 708 WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation 709 enddo !DO iq=1,nqtot 710 WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc 711 WRITE(lunout,*) 'iqfils=',iqfils 712 WRITE(lunout,*) 'nqdesc_tot=',nqdesc_tot 713 714 ! Interdire autres schémas que 10 pour les traceurs fils, et autres schémas 715 ! que 10 et 14 si des pères ont des fils 716 do iq=1,nqtot 717 if (iqpere(iq).gt.0) then 718 ! ce traceur a un père qui n'est pas l'air 719 ! Seul le schéma 10 est autorisé 720 if (iadv(iq)/=10) then 721 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons' 722 CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1) 723 endif 724 ! Le traceur père ne peut être advecté que par schéma 10 ou 14: 725 IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN 726 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for fathers' 727 CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1) 728 endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN 729 endif !if (iqpere(iq).gt.0) the 730 enddo !do iq=1,nqtot 731 732 WRITE(lunout,*) 'infotrac init fin' 733 734 ! detecter quels sont les traceurs isotopiques parmi des traceurs 735 call infotrac_isoinit(tnom_0,nqtrue) 736 737 !----------------------------------------------------------------------- 738 ! Finalize : 739 ! 740 DEALLOCATE(tnom_0, hadv, vadv,tnom_transp) 741 742 743 END SUBROUTINE infotrac_init 744 745 SUBROUTINE infotrac_isoinit(tnom_0,nqtrue) 746 747 #ifdef CPP_IOIPSL 748 use IOIPSL 749 #else 750 ! if not using IOIPSL, we still need to use (a local version of) getin 751 use ioipsl_getincom 752 #endif 753 implicit none 754 755 ! inputs 756 INTEGER nqtrue 757 CHARACTER(len=15) tnom_0(nqtrue) 758 759 ! locals 760 CHARACTER(len=3), DIMENSION(niso_possibles) :: tnom_iso 761 INTEGER, ALLOCATABLE,DIMENSION(:,:) :: nb_iso,nb_traciso 762 INTEGER, ALLOCATABLE,DIMENSION(:) :: nb_isoind 763 INTEGER :: ntraceurs_zone_prec,iq,phase,ixt,iiso,izone 764 CHARACTER(len=19) :: tnom_trac 765 INCLUDE "iniprint.h" 766 767 tnom_iso=(/'eau','HDO','O18','O17','HTO'/) 768 769 ALLOCATE(nb_iso(niso_possibles,nqo)) 770 ALLOCATE(nb_isoind(nqo)) 771 ALLOCATE(nb_traciso(niso_possibles,nqo)) 772 ALLOCATE(iso_num(nqtot)) 773 ALLOCATE(iso_indnum(nqtot)) 774 ALLOCATE(zone_num(nqtot)) 775 ALLOCATE(phase_num(nqtot)) 776 777 iso_num(:)=0 778 iso_indnum(:)=0 779 zone_num(:)=0 780 phase_num(:)=0 781 indnum_fn_num(:)=0 782 use_iso(:)=.false. 783 nb_iso(:,:)=0 784 nb_isoind(:)=0 785 nb_traciso(:,:)=0 786 niso=0 787 ntraceurs_zone=0 788 ntraceurs_zone_prec=0 789 ntraciso=0 790 791 do iq=nqo+1,nqtot 792 ! write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq) 793 do phase=1,nqo 794 do ixt= 1,niso_possibles 795 tnom_trac=trim(tnom_0(phase))//'_' 796 tnom_trac=trim(tnom_trac)//trim(tnom_iso(ixt)) 797 ! write(*,*) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac 798 IF (tnom_0(iq) == tnom_trac) then 799 ! write(lunout,*) 'Ce traceur est un isotope' 800 nb_iso(ixt,phase)=nb_iso(ixt,phase)+1 801 nb_isoind(phase)=nb_isoind(phase)+1 802 iso_num(iq)=ixt 803 iso_indnum(iq)=nb_isoind(phase) 804 indnum_fn_num(ixt)=iso_indnum(iq) 805 phase_num(iq)=phase 806 ! write(lunout,*) 'iso_num(iq)=',iso_num(iq) 807 ! write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq) 808 ! write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt) 809 ! write(lunout,*) 'phase_num(iq)=',phase_num(iq) 810 goto 20 811 else if (iqpere(iq).gt.0) then 812 if (tnom_0(iqpere(iq)) == tnom_trac) then 813 ! write(lunout,*) 'Ce traceur est le fils d''un isotope' 814 ! c'est un traceur d'isotope 815 nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1 816 iso_num(iq)=ixt 817 iso_indnum(iq)=indnum_fn_num(ixt) 818 zone_num(iq)=nb_traciso(ixt,phase) 819 phase_num(iq)=phase 820 ! write(lunout,*) 'iso_num(iq)=',iso_num(iq) 821 ! write(lunout,*) 'phase_num(iq)=',phase_num(iq) 822 ! write(lunout,*) 'zone_num(iq)=',zone_num(iq) 823 goto 20 824 endif !if (tnom_0(iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then 825 endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then 826 enddo !do ixt= niso_possibles 827 enddo !do phase=1,nqo 828 20 continue 829 enddo !do iq=1,nqtot 830 831 ! write(lunout,*) 'iso_num=',iso_num 832 ! write(lunout,*) 'iso_indnum=',iso_indnum 833 ! write(lunout,*) 'zone_num=',zone_num 834 ! write(lunout,*) 'phase_num=',phase_num 835 ! write(lunout,*) 'indnum_fn_num=',indnum_fn_num 836 837 do ixt= 1,niso_possibles 838 839 if (nb_iso(ixt,1).eq.1) then 840 ! on vérifie que toutes les phases ont le même nombre de 841 ! traceurs 842 do phase=2,nqo 843 if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then 844 ! write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase) 845 CALL abort_gcm('infotrac_init','Phases must have same number of isotopes',1) 846 endif 847 enddo !do phase=2,nqo 848 849 niso=niso+1 850 use_iso(ixt)=.true. 851 ntraceurs_zone=nb_traciso(ixt,1) 852 853 ! on vérifie que toutes les phases ont le même nombre de 854 ! traceurs 855 do phase=2,nqo 856 if (nb_traciso(ixt,phase).ne.ntraceurs_zone) then 857 write(lunout,*) 'ixt,phase,nb_traciso=',ixt,phase,nb_traciso(ixt,phase) 858 write(lunout,*) 'ntraceurs_zone=',ntraceurs_zone 859 CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1) 860 endif 861 enddo !do phase=2,nqo 862 ! on vérifie que tous les isotopes ont le même nombre de 863 ! traceurs 864 if (ntraceurs_zone_prec.gt.0) then 865 if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 866 ntraceurs_zone_prec=ntraceurs_zone 867 else !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 868 write(*,*) 'ntraceurs_zone_prec,ntraceurs_zone=',ntraceurs_zone_prec,ntraceurs_zone 869 CALL abort_gcm('infotrac_init', & 870 &'Isotope tracers are not well defined in traceur.def',1) 871 endif !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 872 endif !if (ntraceurs_zone_prec.gt.0) then 873 874 else if (nb_iso(ixt,1).ne.0) then 875 WRITE(lunout,*) 'nqo,ixt=',nqo,ixt 876 WRITE(lunout,*) 'nb_iso(ixt,1)=',nb_iso(ixt,1) 877 CALL abort_gcm('infotrac_init','Isotopes are not well defined in traceur.def',1) 878 endif !if (nb_iso(ixt,1).eq.1) then 879 enddo ! do ixt= niso_possibles 880 881 ! dimensions isotopique: 882 ntraciso=niso*(ntraceurs_zone+1) 883 ! WRITE(lunout,*) 'niso=',niso 884 ! WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso 885 886 ! flags isotopiques: 887 if (niso.gt.0) then 888 ok_isotopes=.true. 889 else 890 ok_isotopes=.false. 891 endif 892 ! WRITE(lunout,*) 'ok_isotopes=',ok_isotopes 893 894 if (ok_isotopes) then 895 ok_iso_verif=.false. 896 call getin('ok_iso_verif',ok_iso_verif) 897 ok_init_iso=.false. 898 call getin('ok_init_iso',ok_init_iso) 899 tnat=(/1.0,155.76e-6,2005.2e-6,0.004/100.,0.0/) 900 alpha_ideal=(/1.0,1.01,1.006,1.003,1.0/) 901 endif !if (ok_isotopes) then 902 ! WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif 903 ! WRITE(lunout,*) 'ok_init_iso=',ok_init_iso 904 905 if (ntraceurs_zone.gt.0) then 906 ok_isotrac=.true. 907 else 908 ok_isotrac=.false. 909 endif 910 ! WRITE(lunout,*) 'ok_isotrac=',ok_isotrac 911 912 ! remplissage du tableau iqiso(ntraciso,phase) 913 ALLOCATE(iqiso(ntraciso,nqo)) 914 iqiso(:,:)=0 915 do iq=1,nqtot 916 if (iso_num(iq).gt.0) then 917 ixt=iso_indnum(iq)+zone_num(iq)*niso 918 iqiso(ixt,phase_num(iq))=iq 919 endif 920 enddo 921 ! WRITE(lunout,*) 'iqiso=',iqiso 922 923 ! replissage du tableau index_trac(ntraceurs_zone,niso) 924 ALLOCATE(index_trac(ntraceurs_zone,niso)) 925 if (ok_isotrac) then 926 do iiso=1,niso 927 do izone=1,ntraceurs_zone 928 index_trac(izone,iiso)=iiso+izone*niso 929 enddo 930 enddo 931 else !if (ok_isotrac) then 932 index_trac(:,:)=0.0 933 endif !if (ok_isotrac) then 934 ! write(lunout,*) 'index_trac=',index_trac 935 936 ! Finalize : 937 DEALLOCATE(nb_iso) 938 939 END SUBROUTINE infotrac_isoinit 528 529 !=== Table used to get iq (index in dyn array, size nqtot) from the isotope and phase indexes ; the full isotopes list 530 ! (including tagging tracers) is sorted this way: iso1, iso2, ..., iso1_zone1, iso2_zone1, ..., iso1_zoneN, iso2_zoneN 531 s%iTraPha = RESHAPE( [( (strIdx(t(:)%name, addPhase(s%trac(it),s%phas(ip:ip))), it=1, s%nitr), ip=1, s%npha)], & 532 [s%nitr, s%npha] ) 533 534 !=== Table used to get ix (index in tagging tracers isotopes list, size nitr) from the zone and isotope indexes 535 s%iZonIso = RESHAPE( [( (strIdx(s%trac(:), TRIM(s%trac(it))//'_'//TRIM(s%zone(iz))), iz=1, s%nzon), it=1, s%niso)], & 536 [s%nzon, s%niso] ) 537 END DO 538 539 !=== Indexes, in dynamical tracers list, of the tracers transmitted to phytrac (nqtottr non-vanishing elements) 540 ll = delPhase(t%name)/='H2O' .AND. t%iso_num ==0 !--- Mask of tracers passed to the physics 541 t(:)%itr = UNPACK([(iq,iq=1,COUNT(ll))], ll, [(0, iq=1, nqtot)]) 542 itr_indice = PACK(t(:)%itr, MASK = t(:)%itr/=0) !--- Might be removed (t%itr should be enough) 543 544 !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE 545 ! DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat AND alpha_ideal) 546 IF(readIsotopesFile('isotopes_params.def',isotopes)) CALL abort_gcm(modname,'Problem when reading isotopes parameters',1) 547 print*,'coincoin' 548 549 !=== Specific to water 550 CALL getKey_init(tracers, isotopes) 551 IF(isoSelect('H2O')) RETURN !--- Select water isotopes ; finished if no water isotopes. 552 iH2O = ixIso !--- Keep track of water family index 553 lerr = getKey('tnat' ,tnat, isoName) 554 lerr = getKey('alpha',alpha_ideal, isoName) 555 CALL msg('end') 556 557 END SUBROUTINE infotrac_init 558 559 560 !============================================================================================================================== 561 !=== THE ROUTINE isoSelect IS USED TO SWITCH FROM AN ISOTOPE FAMILY TO ANOTHER: ISOTOPES DEPENDENT PARAMETERS ARE UPDATED 562 ! Singe generic "isoSelect" routine, using the predefined parent index (fast version) or its name (first time). 563 !============================================================================================================================== 564 LOGICAL FUNCTION isoSelectByName(iName) RESULT(lerr) 565 CHARACTER(LEN=*), INTENT(IN) :: iName 566 INTEGER :: iIso 567 iIso = strIdx(isotopes(:)%prnt, iName) 568 IF(test(fmsg(iIso == 0,'no isotope family named "'//TRIM(iName)//'"'),lerr)) RETURN 569 IF(isoSelectByIndex(iIso)) RETURN 570 END FUNCTION isoSelectByName 571 !============================================================================================================================== 572 LOGICAL FUNCTION isoSelectByIndex(iIso) RESULT(lerr) 573 INTEGER, INTENT(IN) :: iIso 574 lerr = .FALSE. 575 IF(iIso == ixIso) RETURN !--- Nothing to do if the index is already OK 576 IF(test(fmsg(iIso<=0 .OR. iIso>=nbIso,'Inconsistent isotopes family index '//TRIM(int2str(iIso))),lerr)) RETURN 577 ixIso = iIso !--- Update currently selected family index 578 isotope => isotopes(ixIso) !--- Select corresponding component 579 !--- VARIOUS ALIASES 580 isoKeys => isotope%keys; niso = isotope%niso 581 isoName => isotope%trac; nitr = isotope%nitr; isoCheck => isotope%check 582 isoZone => isotope%zone; nzon = isotope%nzon; iZonIso => isotope%iZonIso 583 isoPhas => isotope%phas; npha = isotope%npha; iTraPha => isotope%iTraPha 584 END FUNCTION isoSelectByIndex 585 !============================================================================================================================== 940 586 941 587 END MODULE infotrac -
LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/initdynav.F90
r2622 r3852 6 6 USE IOIPSL 7 7 #endif 8 USE infotrac, ONLY : nqtot , ttext8 USE infotrac, ONLY : nqtot 9 9 use com_io_dyn_mod, only : histaveid,histvaveid,histuaveid, & 10 10 dynhistave_file,dynhistvave_file,dynhistuave_file … … 158 158 159 159 ! DO iq=1,nqtot 160 ! call histdef(histaveid, ttext(iq), ttext(iq), '-', & 160 ! call histdef(histaveid, tracers(iq)%lnam, & 161 ! tracers(iq)%lnam, '-', & 161 162 ! iip1, jjp1, thoriid, llm, 1, llm, zvertiid, & 162 163 ! 32, 'ave(X)', t_ops, t_wrt) -
LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/inithist.F
r2622 r3852 7 7 USE IOIPSL 8 8 #endif 9 USE infotrac, ONLY : nqtot , ttext9 USE infotrac, ONLY : nqtot 10 10 use com_io_dyn_mod, only : histid,histvid,histuid, & 11 11 & dynhist_file,dynhistv_file,dynhistu_file … … 157 157 ! 158 158 ! DO iq=1,nqtot 159 ! call histdef(histid, t text(iq), ttext(iq), '-',160 ! . iip1, jjp1, thoriid, llm, 1, llm, zvertiid,159 ! call histdef(histid, tracers(iq)lnam, tracers(iq)%lnam, 160 ! . '-', iip1, jjp1, thoriid, llm, 1, llm, zvertiid, 161 161 ! . 32, 'inst(X)', t_ops, t_wrt) 162 162 ! enddo -
LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/iso_verif_dyn.F90
r3850 r3852 1 function iso_verif_noNaN_nostop(x,err_msg) 2 implicit none 3 ! si x est NaN, on affiche message 4 ! d'erreur et return 1 si erreur 1 LOGICAL FUNCTION iso_verif_noNaN_nostop(x,err_msg) RESULT(out) 2 USE infotrac, ONLY: isoCheck 3 IMPLICIT NONE 4 !--- Display the message if x is NaN and return .TRUE. if an error occured. 5 REAL, INTENT(IN) :: x 6 CHARACTER(LEN=*), INTENT(IN) :: err_msg 7 include "iniprint.h" 8 REAL, PARAMETER :: borne=1e19 9 out = .FALSE. 10 IF(.NOT.isoCheck) RETURN 11 out = x<=-borne .OR. x>=borne 12 IF(.NOT.out) RETURN 13 WRITE(lunout,*) 'Error detected by iso_verif_noNaN: '//TRIM(err_msg) 14 WRITE(lunout,*) 'x=',x 15 END FUNCTION iso_verif_noNaN_nostop 5 16 6 ! input: 7 real x 8 character*(*) err_msg ! message d''erreur à afficher 17 LOGICAL FUNCTION iso_verif_egalite_nostop(a,b,err_msg) RESULT(out) 18 USE infotrac, ONLY: isoCheck 19 IMPLICIT NONE 20 !--- Display the message if a/=b and return .FALSE. if an error occured. 21 ! Equality is checked for absolute and relative error. 22 REAL, INTENT(IN) :: a,b 23 CHARACTER(LEN=*), INTENT(IN) :: err_msg 24 include "iniprint.h" 25 REAL, PARAMETER :: errmax=1e-8, errmaxrel=1e-3 26 out = .FALSE. 27 IF(.NOT.isoCheck) RETURN 28 out = ABS(a-b)>errmax 29 IF(out) out = ABS((a-b)/MAX(MAX(ABS(b),ABS(a)),1e-18))>errmaxrel 30 IF(.NOT.out) RETURN 31 WRITE(lunout,*) 'Error detected by iso_verif_egalite: '//TRIM(err_msg) 32 WRITE(lunout,*) 'a=',a 33 WRITE(lunout,*) 'b=',b 34 END FUNCTION iso_verif_egalite_nostop 9 35 10 ! output 11 real borne 12 parameter (borne=1e19) 13 integer iso_verif_noNaN_nostop 36 LOGICAL FUNCTION iso_verif_aberrant_nostop(x,kiso,q,err_msg) RESULT(out) 37 USE infotrac, ONLY: isoCheck, tnat 38 IMPLICIT NONE 39 !--- Display the message if a/=b and return .FALSE. if an error occured. 40 ! Equality is checked for absolute and relative error. 41 REAL, INTENT(IN) :: x, q 42 INTEGER, INTENT(IN) :: kiso ! 2=HDO, 1=O18 43 CHARACTER(LEN=*), INTENT(IN) :: err_msg 44 include "iniprint.h" 45 REAL, PARAMETER :: qmin=1e-11, deltaDmax=200.0, deltaDmin=-999.9 46 REAL :: deltaD 47 out = .FALSE. 48 IF(.NOT.isoCheck) RETURN 49 IF(q<qmin) RETURN 50 deltaD = (x/q/tnat(kiso)-1)*1000 51 out = deltaD>deltaDmax .OR. deltaD<deltaDmin 52 IF(.NOT.out) RETURN 53 WRITE(lunout,*) 'Error detected by iso_verif_aberrant: '//TRIM(err_msg) 54 WRITE(lunout,*) 'q = ',q 55 WRITE(lunout,*) 'deltaD = ',deltaD 56 WRITE(lunout,*) 'kiso = ',kiso 57 END FUNCTION iso_verif_aberrant_nostop 14 58 15 if ((x.gt.-borne).and.(x.lt.borne)) then16 iso_verif_noNAN_nostop=017 else18 write(*,*) 'erreur detectee par iso_verif_nonNaN:'19 write(*,*) err_msg20 write(*,*) 'x=',x21 iso_verif_noNaN_nostop=122 endif23 24 return25 end26 27 function iso_verif_egalite_nostop28 : (a,b,err_msg)29 implicit none30 ! compare a et b. Si pas egal, on affiche message31 ! d'erreur et stoppe32 ! pour egalite, on verifie erreur absolue et arreur relative33 34 ! input:35 real a, b36 character*(*) err_msg ! message d''erreur à afficher37 38 ! locals39 real errmax ! erreur maximale en absolu.40 real errmaxrel ! erreur maximale en relatif autorisée41 parameter (errmax=1e-8)42 parameter (errmaxrel=1e-3)43 44 ! output45 integer iso_verif_egalite_nostop46 47 iso_verif_egalite_nostop=048 49 if (abs(a-b).gt.errmax) then50 if (abs((a-b)/max(max(abs(b),abs(a)),1e-18))51 : .gt.errmaxrel) then52 write(*,*) 'erreur detectee par iso_verif_egalite:'53 write(*,*) err_msg54 write(*,*) 'a=',a55 write(*,*) 'b=',b56 iso_verif_egalite_nostop=157 endif58 endif59 60 return61 end62 63 64 function iso_verif_aberrant_nostop65 : (x,iso,q,err_msg)66 USE infotrac67 implicit none68 69 ! input:70 real x,q71 integer iso ! 2=HDO, 1=O1872 character*(*) err_msg ! message d''erreur à afficher73 74 ! locals75 real qmin,deltaD76 real deltaDmax,deltaDmin77 parameter (qmin=1e-11)78 parameter (deltaDmax=200.0,deltaDmin=-999.9)79 80 ! output81 integer iso_verif_aberrant_nostop82 83 iso_verif_aberrant_nostop=084 85 ! verifier que HDO est raisonable86 if (q.gt.qmin) then87 deltaD=(x/q/tnat(iso)-1)*100088 if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then89 write(*,*) 'erreur detectee par iso_verif_aberrant:'90 write(*,*) err_msg91 write(*,*) 'q=',q92 write(*,*) 'deltaD=',deltaD93 write(*,*) 'iso=',iso94 iso_verif_aberrant_nostop=195 endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then96 endif !if (q(i,k,iq).gt.qmin) then97 98 99 return100 end101 -
LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/massbarxy.F90
r2597 r3852 21 21 DO ij=1,ip1jm-1 22 22 massebxy(ij,l)=masse(ij ,l)*alpha2(ij ) + & 23 +masse(ij+1 ,l)*alpha3(ij+1 ) + &24 +masse(ij+iip1,l)*alpha1(ij+iip1) + &25 +masse(ij+iip2,l)*alpha4(ij+iip2)23 masse(ij+1 ,l)*alpha3(ij+1 ) + & 24 masse(ij+iip1,l)*alpha1(ij+iip1) + & 25 masse(ij+iip2,l)*alpha4(ij+iip2) 26 26 END DO 27 27 DO ij=iip1,ip1jm,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO -
LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/writedynav.F90
r2622 r3852 6 6 USE ioipsl 7 7 #endif 8 USE infotrac, ONLY : nqtot , ttext8 USE infotrac, ONLY : nqtot 9 9 use com_io_dyn_mod, only : histaveid, histvaveid, histuaveid 10 10 USE comconst_mod, ONLY: cpp … … 106 106 107 107 ! DO iq=1, nqtot 108 ! call histwrite(histaveid, t text(iq), itau_w, q(:, :, iq), &109 ! iip1*jjp1*llm, ndexu)108 ! call histwrite(histaveid, tracers(iq)%lnam, itau_w, & 109 ! q(:, :, iq), iip1*jjp1*llm, ndexu) 110 110 ! enddo 111 111 -
LMDZ6/branches/LMDZ-tracers/libf/dyn3d_common/writehist.F
r2622 r3852 7 7 USE ioipsl 8 8 #endif 9 USE infotrac, ONLY : nqtot , ttext9 USE infotrac, ONLY : nqtot 10 10 use com_io_dyn_mod, only : histid,histvid,histuid 11 11 USE temps_mod, ONLY: itau_dyn … … 100 100 C 101 101 ! DO iq=1,nqtot 102 ! call histwrite(histid, t text(iq), itau_w, q(:,:,iq),102 ! call histwrite(histid, tracers(iq)%lnam, itau_w, q(:,:,iq), 103 103 ! . iip1*jjp1*llm, ndexu) 104 104 ! enddo
Note: See TracChangeset
for help on using the changeset viewer.
