Context Navigation

infotrac.F90 @ 3998

Last change on this file since 3998 was 3998, checked in by Ehouarn Millour, 3 years ago
Cleanup for tracers. Stop if there is no "traceur.def" file to read, rather than try to guess number and type of tracers. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:keywords set to `Id`
File size: 36.6 KB

Line
1	! $Id: infotrac.F90 3998 2021-11-02 14:54:30Z emillour $
2	!
3	MODULE infotrac
4
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: on retranche les isotopes des traceurs habituels
15	! On fait un tableaux d'indices des traceurs qui passeront dans phytrac
16	INTEGER, SAVE :: nqtottr
17	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: itr_indice
18
19	! CRisi: nb traceurs peres= directement advectes par l'air
20	INTEGER, SAVE :: nqperes
21
22	! ThL: nb traceurs INCA
23	INTEGER, SAVE :: nqINCA
24
25	! ThL: nb traceurs CO2
26	INTEGER, SAVE :: nqCO2
27
28	! Name variables
29	CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
30	CHARACTER(len=23), ALLOCATABLE, DIMENSION(:), SAVE :: ttext ! tracer long name for diagnostics
31
32	! iadv : index of trasport schema for each tracer
33	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iadv
34
35	! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the
36	! dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code.
37	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique
38
39	! CRisi: tableaux de fils
40	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqfils
41	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les generations
42	INTEGER, SAVE :: nqdesc_tot
43	INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqfils
44	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iqpere
45	REAL :: qperemin,masseqmin,ratiomin ! MVals et CRisi
46	PARAMETER (qperemin=1e-16,masseqmin=1e-16,ratiomin=1e-16) ! MVals
47
48	! conv_flg(it)=0 : convection desactivated for tracer number it
49	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: conv_flg
50	! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it
51	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: pbl_flg
52
53	CHARACTER(len=4),SAVE :: type_trac
54	CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym
55
56	! CRisi: cas particulier des isotopes
57	LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso
58	INTEGER :: niso_possibles
59	PARAMETER ( niso_possibles=5)
60	REAL, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal
61	LOGICAL, DIMENSION(niso_possibles),SAVE :: use_iso
62	INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase)
63	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numero iso entre 1 et niso_possibles en fn de nqtot
64	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numero iso entre 1 et niso effectif en fn de nqtot
65	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numero de la zone de tracage en fn de nqtot
66	INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numero de la zone de tracage en fn de nqtot
67	INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numero d isotope entre 1 et niso_possibles
68	INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: index_trac ! numero ixt en fn izone, indnum entre 1 et niso
69	INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso
70
71	#ifdef CPP_StratAer
72	!--CK/OB for stratospheric aerosols
73	INTEGER, SAVE :: nbtr_bin
74	INTEGER, SAVE :: nbtr_sulgas
75	INTEGER, SAVE :: id_OCS_strat
76	INTEGER, SAVE :: id_SO2_strat
77	INTEGER, SAVE :: id_H2SO4_strat
78	INTEGER, SAVE :: id_BIN01_strat
79	INTEGER, SAVE :: id_TEST_strat
80	#endif
81
82	CONTAINS
83
84	SUBROUTINE infotrac_init
85	USE control_mod, ONLY: planet_type, config_inca
86	#ifdef REPROBUS
87	USE CHEM_REP, ONLY : Init_chem_rep_trac
88	#endif
89	IMPLICIT NONE
90	!=======================================================================
91	!
92	! Auteur: P. Le Van /L. Fairhead/F.Hourdin
93	! -------
94	! Modif special traceur F.Forget 05/94
95	! Modif M-A Filiberti 02/02 lecture de traceur.def
96	!
97	! Objet:
98	! ------
99	! GCM LMD nouvelle grille
100	!
101	!=======================================================================
102	! ... modification de l'integration de q ( 26/04/94 ) ....
103	!-----------------------------------------------------------------------
104	! Declarations
105
106	INCLUDE "dimensions.h"
107	INCLUDE "iniprint.h"
108
109	! Local variables
110	INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv ! index of horizontal trasport schema
111	INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv ! index of vertical trasport schema
112
113	INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca ! index of horizontal trasport schema
114	INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca ! index of vertical trasport schema
115
116	INTEGER, ALLOCATABLE, DIMENSION(:) :: conv_flg_inca
117	INTEGER, ALLOCATABLE, DIMENSION(:) :: pbl_flg_inca
118	CHARACTER(len=8), ALLOCATABLE, DIMENSION(:) :: solsym_inca
119
120	CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_0 ! tracer short name
121	CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi
122	CHARACTER(len=3), DIMENSION(30) :: descrq
123	CHARACTER(len=1), DIMENSION(3) :: txts
124	CHARACTER(len=2), DIMENSION(9) :: txtp
125	CHARACTER(len=23) :: str1,str2
126
127	INTEGER :: nqtrue ! number of tracers read from tracer.def, without higer order of moment
128	INTEGER :: iq, new_iq, iiq, jq, ierr,itr
129	INTEGER :: ifils,ipere,generation ! CRisi
130	LOGICAL :: continu,nouveau_traceurdef
131	INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def
132	CHARACTER(len=30) :: tchaine
133
134	character(len=*),parameter :: modname="infotrac_init"
135
136	!-----------------------------------------------------------------------
137	! Initialization :
138	!
139	txts=(/'x','y','z'/)
140	txtp=(/'x ','y ','z ','xx','xy','xz','yy','yz','zz'/)
141
142	descrq(14)='VLH'
143	descrq(10)='VL1'
144	descrq(11)='VLP'
145	descrq(12)='FH1'
146	descrq(13)='FH2'
147	descrq(16)='PPM'
148	descrq(17)='PPS'
149	descrq(18)='PPP'
150	descrq(20)='SLP'
151	descrq(30)='PRA'
152
153
154	! Coherence test between parameter type_trac, config_inca and preprocessing keys
155	IF (type_trac=='inca') THEN
156	WRITE(lunout,*) 'You have chosen to couple with INCA chemistry model : type_trac=', &
157	type_trac,' config_inca=',config_inca
158	IF (config_inca/='aero' .AND. config_inca/='aeNP' .AND. config_inca/='chem') THEN
159	WRITE(lunout,*) 'Incoherence between type_trac and config_inca. Model stops. Modify run.def'
160	CALL abort_gcm('infotrac_init','Incoherence between type_trac and config_inca',1)
161	ENDIF
162	#ifndef INCA
163	WRITE(lunout,*) 'To run this option you must add cpp key INCA and compile with INCA code'
164	CALL abort_gcm('infotrac_init','You must compile with cpp key INCA',1)
165	#endif
166	ELSE IF (type_trac=='repr') THEN
167	WRITE(lunout,*) 'You have chosen to couple with REPROBUS chemestry model : type_trac=', type_trac
168	#ifndef REPROBUS
169	WRITE(lunout,*) 'To run this option you must add cpp key REPROBUS and compile with REPRPBUS code'
170	CALL abort_gcm('infotrac_init','You must compile with cpp key REPROBUS',1)
171	#endif
172	ELSE IF (type_trac == 'co2i') THEN
173	WRITE(lunout,*) 'You have chosen to run with CO2 cycle: type_trac=', type_trac
174	ELSE IF (type_trac == 'coag') THEN
175	WRITE(lunout,*) 'Tracers are treated for COAGULATION tests : type_trac=', type_trac
176	#ifndef CPP_StratAer
177	WRITE(lunout,*) 'To run this option you must add cpp key StratAer and compile with StratAer code'
178	CALL abort_gcm('infotrac_init','You must compile with cpp key StratAer',1)
179	#endif
180	ELSE IF (type_trac == 'lmdz') THEN
181	WRITE(lunout,*) 'Tracers are treated in LMDZ only : type_trac=', type_trac
182	ELSE IF (type_trac == 'inco') THEN ! ThL
183	WRITE(lunout,*) 'Using jointly INCA and CO2 cycle: type_trac =', type_trac
184	IF (config_inca/='aero' .AND. config_inca/='aeNP' .AND. config_inca/='chem') THEN
185	WRITE(lunout,*) 'Incoherence between type_trac and config_inca. Model stops. Modify run.def'
186	CALL abort_gcm('infotrac_init','Incoherence between type_trac and config_inca',1)
187	ENDIF
188	#ifndef INCA
189	WRITE(lunout,*) 'To run this option you must add cpp key INCA and compilewith INCA code'
190	CALL abort_gcm('infotrac_init','You must compile with cpp key INCA',1)
191	#endif
192	ELSE
193	WRITE(lunout,*) 'type_trac=',type_trac,' not possible. Model stops'
194	CALL abort_gcm('infotrac_init','bad parameter',1)
195	ENDIF
196
197	! Test if config_inca is other then none for run without INCA
198	IF (type_trac/='inca' .AND. type_trac/='inco' .AND. config_inca/='none') THEN
199	WRITE(lunout,*) 'config_inca will now be changed to none as you do not couple with INCA model'
200	config_inca='none'
201	ENDIF
202
203	!-----------------------------------------------------------------------
204	!
205	! 1) Get the true number of tracers + water vapor/liquid
206	! Here true tracers (nqtrue) means declared tracers (only first order)
207	!
208	!-----------------------------------------------------------------------
209	IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN
210	IF (type_trac=='co2i') THEN
211	nqCO2 = 1
212	ELSE
213	nqCO2 = 0
214	ENDIF
215	OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
216	IF(ierr.EQ.0) THEN
217	WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
218	READ(90,*) nqtrue
219	write(lunout,*) 'nqtrue=',nqtrue
220	ELSE
221	WRITE(lunout,*) trim(modname),': Failed opening traceur.def'
222	CALL abort_gcm(modname,"file traceur.def not found!",1)
223	ENDIF
224	!jyg<
225	!! if ( planet_type=='earth') then
226	!! ! For Earth, water vapour & liquid tracers are not in the physics
227	!! nbtr=nqtrue-2
228	!! else
229	!! ! Other planets (for now); we have the same number of tracers
230	!! ! in the dynamics than in the physics
231	!! nbtr=nqtrue
232	!! endif
233	!>jyg
234	ELSE ! type_trac=inca or inco
235	IF (type_trac=='inco') THEN
236	nqCO2 = 1
237	ELSE
238	nqCO2 = 0
239	ENDIF
240	!jyg<
241	! The traceur.def file is used to define the number "nqo" of water phases
242	! present in the simulation. Default : nqo = 2.
243	OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
244	IF(ierr.EQ.0) THEN
245	WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
246	READ(90,*) nqo
247	ELSE
248	WRITE(lunout,*) trim(modname),': Failed opening traceur.def'
249	CALL abort_gcm(modname,"file traceur.def not found!",1)
250	ENDIF
251	IF (nqo /= 2 .AND. nqo /= 3 ) THEN
252	IF (nqo == 4 .AND. type_trac=='inco') THEN ! ThL
253	WRITE(lunout,*) trim(modname),': you are coupling with INCA, and also using CO2i.'
254	nqo = 3 ! A ameliorier... je force 3 traceurs eau... ThL
255	WRITE(lunout,*) trim(modname),': nqo = ',nqo
256	ELSE
257	WRITE(lunout,*) trim(modname),': nqo=',nqo, ' is not allowed. Only 2 or 3 water phases allowed'
258	CALL abort_gcm('infotrac_init','Bad number of water phases',1)
259	ENDIF
260	ENDIF
261	! nbtr has been read from INCA by init_const_lmdz() in gcm.F
262	#ifdef INCA
263	CALL Init_chem_inca_trac(nqINCA)
264	#else
265	nqINCA=0
266	#endif
267	nbtr=nqINCA+nqCO2
268	nqtrue=nbtr+nqo
269	WRITE(lunout,*) trim(modname),': nqo = ',nqo
270	WRITE(lunout,*) trim(modname),': nbtr = ',nbtr
271	WRITE(lunout,*) trim(modname),': nqtrue = ',nqtrue
272	WRITE(lunout,*) trim(modname),': nqCO2 = ',nqCO2
273	WRITE(lunout,*) trim(modname),': nqINCA = ',nqINCA
274	ALLOCATE(hadv_inca(nqINCA), vadv_inca(nqINCA), conv_flg_inca(nqINCA), pbl_flg_inca(nqINCA), solsym_inca(nqINCA))
275	ENDIF ! type_trac 'inca' ou 'inco'
276	!>jyg
277
278	IF ((planet_type=="earth").and.(nqtrue < 2)) THEN
279	WRITE(lunout,*) trim(modname),': nqtrue=',nqtrue, ' is not allowed. 2 tracers is the minimum'
280	CALL abort_gcm('infotrac_init','Not enough tracers',1)
281	ENDIF
282
283	!jyg<
284
285	!
286	! Allocate variables depending on nqtrue
287	!
288	ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue))
289
290
291	!-----------------------------------------------------------------------
292	! 2) Choix des schemas d'advection pour l'eau et les traceurs
293	!
294	! iadv = 1 schema transport type "humidite specifique LMD"
295	! iadv = 2 schema amont
296	! iadv = 14 schema Van-leer + humidite specifique
297	! Modif F.Codron
298	! iadv = 10 schema Van-leer (retenu pour l'eau vapeur et liquide)
299	! iadv = 11 schema Van-Leer pour hadv et version PPM (Monotone) pour vadv
300	! iadv = 12 schema Frederic Hourdin I
301	! iadv = 13 schema Frederic Hourdin II
302	! iadv = 16 schema PPM Monotone(Collela & Woodward 1984)
303	! iadv = 17 schema PPM Semi Monotone (overshoots autorises)
304	! iadv = 18 schema PPM Positif Defini (overshoots undershoots autorises)
305	! iadv = 20 schema Slopes
306	! iadv = 30 schema Prather
307	!
308	! Dans le tableau q(ij,l,iq) : iq = 1 pour l'eau vapeur
309	! iq = 2 pour l'eau liquide
310	! Et eventuellement iq = 3,nqtot pour les autres traceurs
311	!
312	! iadv(1): choix pour l'eau vap. et iadv(2) : choix pour l'eau liq.
313	!------------------------------------------------------------------------
314	!
315	! Get choice of advection schema from file tracer.def or from INCA
316	!---------------------------------------------------------------------
317	IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN
318
319	! Continue to read tracer.def
320	DO iq=1,nqtrue
321
322	write(,) 'infotrac 237: iq=',iq
323	! CRisi: ajout du nom du fluide transporteur
324	! mais rester retro compatible
325	READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
326	write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
327	write(lunout,*) 'tchaine=',trim(tchaine)
328	write(,) 'infotrac 238: IOstatus=',IOstatus
329	if (IOstatus.ne.0) then
330	CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
331	endif
332	! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
333	! espace ou pas au milieu de la chaine.
334	continu=.true.
335	nouveau_traceurdef=.false.
336	iiq=1
337	do while (continu)
338	if (tchaine(iiq:iiq).eq.' ') then
339	nouveau_traceurdef=.true.
340	continu=.false.
341	else if (iiq.lt.LEN_TRIM(tchaine)) then
342	iiq=iiq+1
343	else
344	continu=.false.
345	endif
346	enddo
347	write(,) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
348	if (nouveau_traceurdef) then
349	write(lunout,*) 'C''est la nouvelle version de traceur.def'
350	tnom_0(iq)=tchaine(1:iiq-1)
351	tnom_transp(iq)=tchaine(iiq+1:30)
352	else
353	write(lunout,*) 'C''est l''ancienne version de traceur.def'
354	write(lunout,*) 'On suppose que les traceurs sont tous d''air'
355	tnom_0(iq)=tchaine
356	tnom_transp(iq) = 'air'
357	endif
358	write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
359	write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
360
361	ENDDO!DO iq=1,nqtrue
362
363	CLOSE(90)
364
365	WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
366	WRITE(lunout,*) trim(modname),': nombre total de traceurs ',nqtrue
367	DO iq=1,nqtrue
368	WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq)
369	END DO
370
371	IF ( planet_type=='earth') THEN
372	!CR: nombre de traceurs de l eau
373	IF (tnom_0(3) == 'H2Oi') THEN
374	nqo=3
375	ELSE
376	nqo=2
377	ENDIF
378	! For Earth, water vapour & liquid tracers are not in the physics
379	nbtr=nqtrue-nqo
380	ELSE
381	! Other planets (for now); we have the same number of tracers
382	! in the dynamics than in the physics
383	nbtr=nqtrue
384	ENDIF
385
386	#ifdef CPP_StratAer
387	IF (type_trac == 'coag') THEN
388	nbtr_bin=0
389	nbtr_sulgas=0
390	DO iq=1,nqtrue
391	IF (tnom_0(iq)(1:3)=='BIN') THEN !check if tracer name contains 'BIN'
392	nbtr_bin=nbtr_bin+1
393	ENDIF
394	IF (tnom_0(iq)(1:3)=='GAS') THEN !check if tracer name contains 'GAS'
395	nbtr_sulgas=nbtr_sulgas+1
396	ENDIF
397	ENDDO
398	print*,'nbtr_bin=',nbtr_bin
399	print*,'nbtr_sulgas=',nbtr_sulgas
400	DO iq=1,nqtrue
401	IF (tnom_0(iq)=='GASOCS') THEN
402	id_OCS_strat=iq-nqo
403	ENDIF
404	IF (tnom_0(iq)=='GASSO2') THEN
405	id_SO2_strat=iq-nqo
406	ENDIF
407	IF (tnom_0(iq)=='GASH2SO4') THEN
408	id_H2SO4_strat=iq-nqo
409	ENDIF
410	IF (tnom_0(iq)=='BIN01') THEN
411	id_BIN01_strat=iq-nqo
412	ENDIF
413	IF (tnom_0(iq)=='GASTEST') THEN
414	id_TEST_strat=iq-nqo
415	ENDIF
416	ENDDO
417	print*,'id_OCS_strat =',id_OCS_strat
418	print*,'id_SO2_strat =',id_SO2_strat
419	print*,'id_H2SO4_strat=',id_H2SO4_strat
420	print*,'id_BIN01_strat=',id_BIN01_strat
421	ENDIF
422	#endif
423
424	ENDIF ! (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac = 'coag' .OR. type_trac = 'co2i')
425	!jyg<
426	!
427
428	! Transfert number of tracers to Reprobus
429	IF (type_trac == 'repr') THEN
430	#ifdef REPROBUS
431	CALL Init_chem_rep_trac(nbtr,nqo,tnom_0)
432	#endif
433	ENDIF
434	!
435	! Allocate variables depending on nbtr
436	!
437	ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
438	conv_flg(:) = 1 ! convection activated for all tracers
439	pbl_flg(:) = 1 ! boundary layer activated for all tracers
440
441	IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN ! config_inca='aero' ou 'chem'
442	!>jyg
443	! le module de chimie fournit les noms des traceurs
444	! et les schemas d'advection associes. excepte pour ceux lus
445	! dans traceur.def
446
447	DO iq=1,nqo+nqCO2
448
449	write(,) 'infotrac 237: iq=',iq
450	! CRisi: ajout du nom du fluide transporteur
451	! mais rester retro compatible
452	READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
453	write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
454	write(lunout,*) 'tchaine=',trim(tchaine)
455	write(,) 'infotrac 238: IOstatus=',IOstatus
456	if (IOstatus.ne.0) then
457	CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
458	endif
459	! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
460	! espace ou pas au milieu de la chaine.
461	continu=.true.
462	nouveau_traceurdef=.false.
463	iiq=1
464
465	do while (continu)
466	if (tchaine(iiq:iiq).eq.' ') then
467	nouveau_traceurdef=.true.
468	continu=.false.
469	else if (iiq.lt.LEN_TRIM(tchaine)) then
470	iiq=iiq+1
471	else
472	continu=.false.
473	endif
474	enddo
475
476	write(,) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
477
478	if (nouveau_traceurdef) then
479	write(lunout,*) 'C''est la nouvelle version de traceur.def'
480	tnom_0(iq)=tchaine(1:iiq-1)
481	tnom_transp(iq)=tchaine(iiq+1:30)
482	else
483	write(lunout,*) 'C''est l''ancienne version de traceur.def'
484	write(lunout,*) 'On suppose que les traceurs sont tous d''air'
485	tnom_0(iq)=tchaine
486	tnom_transp(iq) = 'air'
487	endif
488
489	write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
490	write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
491
492	ENDDO !DO iq=1,nqo
493	CLOSE(90)
494
495
496	#ifdef INCA
497	CALL init_transport( &
498	hadv_inca, &
499	vadv_inca, &
500	conv_flg_inca, &
501	pbl_flg_inca, &
502	solsym_inca)
503
504	conv_flg(1+nqCO2:nbtr) = conv_flg_inca
505	pbl_flg(1+nqCO2:nbtr) = pbl_flg_inca
506	solsym(1+nqCO2:nbtr) = solsym_inca
507
508	IF (type_trac == 'inco') THEN
509	conv_flg(1:nqCO2) = 1
510	pbl_flg(1:nqCO2) = 1
511	solsym(1:nqCO2) = 'CO2'
512	ENDIF
513	#endif
514
515	!jyg<
516	DO iq = nqo+nqCO2+1, nqtrue
517	hadv(iq) = hadv_inca(iq-nqo-nqCO2)
518	vadv(iq) = vadv_inca(iq-nqo-nqCO2)
519	tnom_0(iq)=solsym_inca(iq-nqo-nqCO2)
520	tnom_transp(iq) = 'air'
521	END DO
522
523	ENDIF ! (type_trac == 'inca' or 'inco')
524
525	!-----------------------------------------------------------------------
526	!
527	! 3) Verify if advection schema 20 or 30 choosen
528	! Calculate total number of tracers needed: nqtot
529	! Allocate variables depending on total number of tracers
530	!-----------------------------------------------------------------------
531	new_iq=0
532	DO iq=1,nqtrue
533	! Add tracers for certain advection schema
534	IF (hadv(iq)<20 .AND. vadv(iq)<20 ) THEN
535	new_iq=new_iq+1 ! no tracers added
536	ELSE IF (hadv(iq)==20 .AND. vadv(iq)==20 ) THEN
537	new_iq=new_iq+4 ! 3 tracers added
538	ELSE IF (hadv(iq)==30 .AND. vadv(iq)==30 ) THEN
539	new_iq=new_iq+10 ! 9 tracers added
540	ELSE
541	WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
542	CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
543	ENDIF
544	END DO
545
546	IF (new_iq /= nqtrue) THEN
547	! The choice of advection schema imposes more tracers
548	! Assigne total number of tracers
549	nqtot = new_iq
550
551	WRITE(lunout,*) trim(modname),': The choice of advection schema for one or more tracers'
552	WRITE(lunout,*) 'makes it necessary to add tracers'
553	WRITE(lunout,*) trim(modname)//': ',nqtrue,' is the number of true tracers'
554	WRITE(lunout,*) trim(modname)//': ',nqtot, ' is the total number of tracers needed'
555
556	ELSE
557	! The true number of tracers is also the total number
558	nqtot = nqtrue
559	ENDIF
560
561	!
562	! Allocate variables with total number of tracers, nqtot
563	!
564	ALLOCATE(tname(nqtot), ttext(nqtot))
565	ALLOCATE(iadv(nqtot), niadv(nqtot))
566
567	!-----------------------------------------------------------------------
568	!
569	! 4) Determine iadv, long and short name
570	!
571	!-----------------------------------------------------------------------
572	new_iq=0
573	DO iq=1,nqtrue
574	new_iq=new_iq+1
575
576	! Verify choice of advection schema
577	IF (hadv(iq)==vadv(iq)) THEN
578	iadv(new_iq)=hadv(iq)
579	ELSE IF (hadv(iq)==10 .AND. vadv(iq)==16) THEN
580	iadv(new_iq)=11
581	ELSE
582	WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
583
584	CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
585	ENDIF
586
587	str1=tnom_0(iq)
588	tname(new_iq)= tnom_0(iq)
589	IF (iadv(new_iq)==0) THEN
590	ttext(new_iq)=trim(str1)
591	ELSE
592	ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
593	ENDIF
594
595	! schemas tenant compte des moments d'ordre superieur
596	str2=ttext(new_iq)
597	IF (iadv(new_iq)==20) THEN
598	DO jq=1,3
599	new_iq=new_iq+1
600	iadv(new_iq)=-20
601	ttext(new_iq)=trim(str2)//txts(jq)
602	tname(new_iq)=trim(str1)//txts(jq)
603	END DO
604	ELSE IF (iadv(new_iq)==30) THEN
605	DO jq=1,9
606	new_iq=new_iq+1
607	iadv(new_iq)=-30
608	ttext(new_iq)=trim(str2)//txtp(jq)
609	tname(new_iq)=trim(str1)//txtp(jq)
610	END DO
611	ENDIF
612	END DO
613
614	!
615	! Find vector keeping the correspodence between true and total tracers
616	!
617	niadv(:)=0
618	iiq=0
619	DO iq=1,nqtot
620	IF(iadv(iq).GE.0) THEN
621	! True tracer
622	iiq=iiq+1
623	niadv(iiq)=iq
624	ENDIF
625	END DO
626
627
628	WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
629	WRITE(lunout,*) trim(modname),': iadv niadv tname ttext :'
630
631	DO iq=1,nqtot
632	WRITE(lunout,*) iadv(iq),niadv(iq), ' ',trim(tname(iq)),' ',trim(ttext(iq))
633	END DO
634
635	!
636	! Test for advection schema.
637	! This version of LMDZ only garantees iadv=10 and iadv=14 (14 only for water vapour) .
638	!
639	DO iq=1,nqtot
640	IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN
641	WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
642	CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1)
643	ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN
644	WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
645	CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
646	ENDIF
647	END DO
648
649
650	! CRisi: quels sont les traceurs fils et les traceurs peres.
651	! initialiser tous les tableaux d'indices lies aux traceurs familiaux
652	! + verifier que tous les peres sont ecrits en premieres positions
653	ALLOCATE(nqfils(nqtot),nqdesc(nqtot))
654	ALLOCATE(iqfils(nqtot,nqtot))
655	ALLOCATE(iqpere(nqtot))
656	nqperes=0
657	nqfils(:)=0
658	nqdesc(:)=0
659	iqfils(:,:)=0
660	iqpere(:)=0
661	nqdesc_tot=0
662	DO iq=1,nqtot
663	if (tnom_transp(iq) == 'air') then
664	! ceci est un traceur pere
665	WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere'
666	nqperes=nqperes+1
667	iqpere(iq)=0
668	else !if (tnom_transp(iq) == 'air') then
669	! ceci est un fils. Qui est son pere?
670	WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un fils'
671	continu=.true.
672	ipere=1
673	do while (continu)
674	if (tnom_transp(iq) == tnom_0(ipere)) then
675	! Son pere est ipere
676	WRITE(lunout,*) 'Le traceur',iq,'appele ', &
677	& trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere))
678	if (iq.eq.ipere) then
679	CALL abort_gcm('infotrac_init','Un fils est son propre pere',1)
680	endif
681	nqfils(ipere)=nqfils(ipere)+1
682	iqfils(nqfils(ipere),ipere)=iq
683	iqpere(iq)=ipere
684	continu=.false.
685	else !if (tnom_transp(iq) == tnom_0(ipere)) then
686	ipere=ipere+1
687	if (ipere.gt.nqtot) then
688	WRITE(lunout,*) 'Le traceur',iq,'appele ', &
689	& trim(tnom_0(iq)),', est orphelin.'
690	CALL abort_gcm('infotrac_init','Un traceur est orphelin',1)
691	endif !if (ipere.gt.nqtot) then
692	endif !if (tnom_transp(iq) == tnom_0(ipere)) then
693	enddo !do while (continu)
694	endif !if (tnom_transp(iq) == 'air') then
695	enddo !DO iq=1,nqtot
696	WRITE(lunout,*) 'infotrac: nqperes=',nqperes
697	WRITE(lunout,*) 'nqfils=',nqfils
698	WRITE(lunout,*) 'iqpere=',iqpere
699	WRITE(lunout,*) 'iqfils=',iqfils
700
701	! Calculer le nombre de descendants a partir de iqfils et de nbfils
702	DO iq=1,nqtot
703	generation=0
704	continu=.true.
705	ifils=iq
706	do while (continu)
707	ipere=iqpere(ifils)
708	if (ipere.gt.0) then
709	nqdesc(ipere)=nqdesc(ipere)+1
710	nqdesc_tot=nqdesc_tot+1
711	iqfils(nqdesc(ipere),ipere)=iq
712	ifils=ipere
713	generation=generation+1
714	else !if (ipere.gt.0) then
715	continu=.false.
716	endif !if (ipere.gt.0) then
717	enddo !do while (continu)
718	WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation
719	enddo !DO iq=1,nqtot
720	WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc
721	WRITE(lunout,*) 'iqfils=',iqfils
722	WRITE(lunout,*) 'nqdesc_tot=',nqdesc_tot
723
724	! Interdire autres schemas que 10 pour les traceurs fils, et autres schemas
725	! que 10 et 14 si des peres ont des fils
726	do iq=1,nqtot
727	if (iqpere(iq).gt.0) then
728	! ce traceur a un pere qui n'est pas l'air
729	! Seul le schema 10 est autorise
730	if (iadv(iq)/=10) then
731	WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons'
732	CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1)
733	endif
734	! Le traceur pere ne peut etre advecte que par schema 10 ou 14:
735	IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
736	WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for fathers'
737	CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1)
738	endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
739	endif !if (iqpere(iq).gt.0) the
740	enddo !do iq=1,nqtot
741
742
743
744	! detecter quels sont les traceurs isotopiques parmi des traceurs
745	call infotrac_isoinit(tnom_0,nqtrue)
746
747	! if (ntraciso.gt.0) then
748	! le 18 sep 2020: on enleve la condition ntraciso.gt.0 car nqtottr doit etre
749	! connu meme si il n'y a pas d'isotopes!
750	write(lunout,*) 'infotrac 702: nbtr,ntraciso=',nbtr,ntraciso
751	! retrancher les traceurs isotopiques de la liste des traceurs qui passent dans
752	! phytrac
753	nbtr=nbtr-nqo*ntraciso
754
755	! faire un tableau d'indice des traceurs qui passeront dans phytrac
756	nqtottr=nqtot-nqo*(1+ntraciso)
757	write(lunout,*) 'infotrac 704: nqtottr,nqtot,nqo=',nqtottr,nqtot,nqo
758	! Rq: nqtottr n'est pas forcement egal a nbtr dans le cas ou new_iq /= nqtrue
759	ALLOCATE (itr_indice(nqtottr))
760	itr_indice(:)=0
761	itr=0
762	do iq=nqo+1, nqtot
763	if (iso_num(iq).eq.0) then
764	itr=itr+1
765	write(,) 'itr=',itr
766	itr_indice(itr)=iq
767	endif !if (iso_num(iq).eq.0) then
768	enddo
769	if (itr.ne.nqtottr) then
770	CALL abort_gcm('infotrac_init','pb dans le calcul de nqtottr',1)
771	endif
772	write(lunout,*) 'itr_indice=',itr_indice
773	! endif !if (ntraciso.gt.0) then
774
775	!-----------------------------------------------------------------------
776	! Finalize :
777	!
778	DEALLOCATE(tnom_0, hadv, vadv,tnom_transp)
779
780	WRITE(lunout,*) 'infotrac init fin'
781
782	END SUBROUTINE infotrac_init
783
784	SUBROUTINE infotrac_isoinit(tnom_0,nqtrue)
785
786	#ifdef CPP_IOIPSL
787	use IOIPSL
788	#else
789	! if not using IOIPSL, we still need to use (a local version of) getin
790	use ioipsl_getincom
791	#endif
792	implicit none
793
794	! inputs
795	INTEGER nqtrue
796	CHARACTER(len=30) tnom_0(nqtrue)
797
798	! locals
799	CHARACTER(len=3), DIMENSION(niso_possibles) :: tnom_iso
800	INTEGER, ALLOCATABLE,DIMENSION(:,:) :: nb_iso,nb_traciso
801	INTEGER, ALLOCATABLE,DIMENSION(:) :: nb_isoind
802	INTEGER :: ntraceurs_zone_prec,iq,phase,ixt,iiso,izone
803	CHARACTER(len=19) :: tnom_trac
804	INCLUDE "iniprint.h"
805
806	tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
807
808	ALLOCATE(nb_iso(niso_possibles,nqo))
809	ALLOCATE(nb_isoind(nqo))
810	ALLOCATE(nb_traciso(niso_possibles,nqo))
811	ALLOCATE(iso_num(nqtot))
812	ALLOCATE(iso_indnum(nqtot))
813	ALLOCATE(zone_num(nqtot))
814	ALLOCATE(phase_num(nqtot))
815
816	iso_num(:)=0
817	iso_indnum(:)=0
818	zone_num(:)=0
819	phase_num(:)=0
820	indnum_fn_num(:)=0
821	use_iso(:)=.false.
822	nb_iso(:,:)=0
823	nb_isoind(:)=0
824	nb_traciso(:,:)=0
825	niso=0
826	ntraceurs_zone=0
827	ntraceurs_zone_prec=0
828	ntraciso=0
829
830	do iq=nqo+1,nqtot
831	! write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq)
832	do phase=1,nqo
833	do ixt= 1,niso_possibles
834	tnom_trac=trim(tnom_0(phase))//'_'
835	tnom_trac=trim(tnom_trac)//trim(tnom_iso(ixt))
836	! write(,) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac
837	IF (tnom_0(iq) == tnom_trac) then
838	! write(lunout,*) 'Ce traceur est un isotope'
839	nb_iso(ixt,phase)=nb_iso(ixt,phase)+1
840	nb_isoind(phase)=nb_isoind(phase)+1
841	iso_num(iq)=ixt
842	iso_indnum(iq)=nb_isoind(phase)
843	indnum_fn_num(ixt)=iso_indnum(iq)
844	phase_num(iq)=phase
845	! write(lunout,*) 'iso_num(iq)=',iso_num(iq)
846	! write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq)
847	! write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt)
848	! write(lunout,*) 'phase_num(iq)=',phase_num(iq)
849	goto 20
850	else if (iqpere(iq).gt.0) then
851	if (tnom_0(iqpere(iq)) == tnom_trac) then
852	! write(lunout,*) 'Ce traceur est le fils d''un isotope'
853	! c'est un traceur d'isotope
854	nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1
855	iso_num(iq)=ixt
856	iso_indnum(iq)=indnum_fn_num(ixt)
857	zone_num(iq)=nb_traciso(ixt,phase)
858	phase_num(iq)=phase
859	! write(lunout,*) 'iso_num(iq)=',iso_num(iq)
860	! write(lunout,*) 'phase_num(iq)=',phase_num(iq)
861	! write(lunout,*) 'zone_num(iq)=',zone_num(iq)
862	goto 20
863	endif !if (tnom_0(iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
864	endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
865	enddo !do ixt= niso_possibles
866	enddo !do phase=1,nqo
867	20 continue
868	enddo !do iq=1,nqtot
869
870	! write(lunout,*) 'iso_num=',iso_num
871	! write(lunout,*) 'iso_indnum=',iso_indnum
872	! write(lunout,*) 'zone_num=',zone_num
873	! write(lunout,*) 'phase_num=',phase_num
874	! write(lunout,*) 'indnum_fn_num=',indnum_fn_num
875
876	do ixt= 1,niso_possibles
877
878	if (nb_iso(ixt,1).eq.1) then
879	! on verifie que toutes les phases ont le meme nombre de
880	! traceurs
881	do phase=2,nqo
882	if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then
883	! write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase)
884	CALL abort_gcm('infotrac_init','Phases must have same number of isotopes',1)
885	endif
886	enddo !do phase=2,nqo
887
888	niso=niso+1
889	use_iso(ixt)=.true.
890	ntraceurs_zone=nb_traciso(ixt,1)
891
892	! on verifie que toutes les phases ont le meme nombre de
893	! traceurs
894	do phase=2,nqo
895	if (nb_traciso(ixt,phase).ne.ntraceurs_zone) then
896	write(lunout,*) 'ixt,phase,nb_traciso=',ixt,phase,nb_traciso(ixt,phase)
897	write(lunout,*) 'ntraceurs_zone=',ntraceurs_zone
898	CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1)
899	endif
900	enddo !do phase=2,nqo
901	! on verifie que tous les isotopes ont le meme nombre de
902	! traceurs
903	if (ntraceurs_zone_prec.gt.0) then
904	if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
905	ntraceurs_zone_prec=ntraceurs_zone
906	else !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
907	write(,) 'ntraceurs_zone_prec,ntraceurs_zone=',ntraceurs_zone_prec,ntraceurs_zone
908	CALL abort_gcm('infotrac_init', &
909	&'Isotope tracers are not well defined in traceur.def',1)
910	endif !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
911	endif !if (ntraceurs_zone_prec.gt.0) then
912
913	else if (nb_iso(ixt,1).ne.0) then
914	WRITE(lunout,*) 'nqo,ixt=',nqo,ixt
915	WRITE(lunout,*) 'nb_iso(ixt,1)=',nb_iso(ixt,1)
916	CALL abort_gcm('infotrac_init','Isotopes are not well defined in traceur.def',1)
917	endif !if (nb_iso(ixt,1).eq.1) then
918	enddo ! do ixt= niso_possibles
919
920	! dimensions isotopique:
921	ntraciso=niso*(ntraceurs_zone+1)
922	! WRITE(lunout,*) 'niso=',niso
923	! WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso
924
925	! flags isotopiques:
926	if (niso.gt.0) then
927	ok_isotopes=.true.
928	else
929	ok_isotopes=.false.
930	endif
931	! WRITE(lunout,*) 'ok_isotopes=',ok_isotopes
932
933	if (ok_isotopes) then
934	ok_iso_verif=.false.
935	call getin('ok_iso_verif',ok_iso_verif)
936	ok_init_iso=.false.
937	call getin('ok_init_iso',ok_init_iso)
938	tnat=(/1.0,155.76e-6,2005.2e-6,0.004/100.,0.0/)
939	alpha_ideal=(/1.0,1.01,1.006,1.003,1.0/)
940	endif !if (ok_isotopes) then
941	! WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif
942	! WRITE(lunout,*) 'ok_init_iso=',ok_init_iso
943
944	if (ntraceurs_zone.gt.0) then
945	ok_isotrac=.true.
946	else
947	ok_isotrac=.false.
948	endif
949	! WRITE(lunout,*) 'ok_isotrac=',ok_isotrac
950
951	! remplissage du tableau iqiso(ntraciso,phase)
952	ALLOCATE(iqiso(ntraciso,nqo))
953	iqiso(:,:)=0
954	do iq=1,nqtot
955	if (iso_num(iq).gt.0) then
956	ixt=iso_indnum(iq)+zone_num(iq)*niso
957	iqiso(ixt,phase_num(iq))=iq
958	endif
959	enddo
960	! WRITE(lunout,*) 'iqiso=',iqiso
961
962	! replissage du tableau index_trac(ntraceurs_zone,niso)
963	ALLOCATE(index_trac(ntraceurs_zone,niso))
964	if (ok_isotrac) then
965	do iiso=1,niso
966	do izone=1,ntraceurs_zone
967	index_trac(izone,iiso)=iiso+izone*niso
968	enddo
969	enddo
970	else !if (ok_isotrac) then
971	index_trac(:,:)=0.0
972	endif !if (ok_isotrac) then
973	! write(lunout,*) 'index_trac=',index_trac
974
975	! Finalize :
976	DEALLOCATE(nb_iso)
977
978	END SUBROUTINE infotrac_isoinit
979
980	END MODULE infotrac

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format