Context Navigation

infotrac.F90 @ 2300

Last change on this file since 2300 was 2270, checked in by crisi, 10 years ago
Adding isotopes in the dynamics and more generally tracers of tracers. CRisi
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
File size: 30.1 KB

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

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

Download in other formats:

Original Format