Context Navigation

infotrac.F90 @ 3888

Last change on this file since 3888 was 3872, checked in by oboucher, 4 years ago
Last modifications for the CO2/INCA (inco) case
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.9 KB

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

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

Download in other formats:

Original Format