Context Navigation

traclmdz_mod.F90 @ 4056

Last change on this file since 4056 was 4056, checked in by dcugnet, 2 years ago

Most of the changes are intended to help to eventually remove the constraints about the tracers assumptions, in particular water tracers.

Remove index tables itr_indice and niadv, replaced by tracers(:)%isAdvected and tracers(:)%isH2OFamily. Most of the loops are now from 1 to nqtot:
- DO iq=nqo+1,nqtot loops are replaced with: DO iq=1,nqtot

IF(tracers(iq)%isH2Ofamily) CYCLE

DO it=1,nbtr; iq=niadv(it+nqo)

and DO it=1,nqtottr; iq=itr_indice(it) loops are replaced with:

it = 0
DO iq = 1, nqtot

IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
it = it+1

Move some StratAer? related code from infotrac to infotrac_phy
Remove "nqperes" variable:

DO iq=1,nqpere loops are replaced with:
DO iq=1,nqtot

IF(tracers(iq)%parent/='air') CYCLE

Cosmetic changes (justification, SELECT CASE instead of multiple IF...) mostly in advtrac* routines.

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: 22.1 KB

Line
1	!$Id $
2	!
3	MODULE traclmdz_mod
4
5	!
6	! In this module all tracers specific to LMDZ are treated. This module is used
7	! only if running without any other chemestry model as INCA or REPROBUS.
8	!
9	IMPLICIT NONE
10
11	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: masktr ! Masque reservoir de sol traceur
12	!$OMP THREADPRIVATE(masktr) ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 )
13	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: fshtr ! Flux surfacique dans le reservoir de sol
14	!$OMP THREADPRIVATE(fshtr)
15	REAL,DIMENSION(:),ALLOCATABLE,SAVE :: hsoltr ! Epaisseur equivalente du reservoir de sol
16	!$OMP THREADPRIVATE(hsoltr)
17	!
18	!Radioelements:
19	!--------------
20	!
21	REAL,DIMENSION(:),ALLOCATABLE,SAVE :: tautr ! Constante de decroissance radioactive
22	!$OMP THREADPRIVATE(tautr)
23	REAL,DIMENSION(:),ALLOCATABLE,SAVE :: vdeptr ! Vitesse de depot sec dans la couche Brownienne
24	!$OMP THREADPRIVATE(vdeptr)
25	REAL,DIMENSION(:),ALLOCATABLE,SAVE :: scavtr ! Coefficient de lessivage
26	!$OMP THREADPRIVATE(scavtr)
27	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srcbe ! Production du beryllium7 dans l atmosphere (U/s/kgA)
28	!$OMP THREADPRIVATE(srcbe)
29
30	LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: radio ! radio(it) = true => decroisssance radioactive
31	!$OMP THREADPRIVATE(radio)
32
33	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: trs ! Conc. radon ds le sol
34	!$OMP THREADPRIVATE(trs)
35
36	INTEGER,SAVE :: id_aga ! Identification number for tracer : Age of stratospheric air
37	!$OMP THREADPRIVATE(id_aga)
38	INTEGER,SAVE :: lev_1p5km ! Approximative vertical layer number at 1.5km above surface, used for calculation of the age of air. The result shouldn't be that sensible to the exactness of this value as long as it is in the lower troposphere.
39	!$OMP THREADPRIVATE(lev_1p5km)
40
41	INTEGER,SAVE :: id_rn, id_pb ! Identification number for tracer : radon (Rn222), lead (Pb210)
42	!$OMP THREADPRIVATE(id_rn, id_pb)
43
44	INTEGER,SAVE :: id_be ! Activation et position du traceur Be7 [ id_be=0 -> desactive ]
45	!$OMP THREADPRIVATE(id_be)
46
47	INTEGER,SAVE :: id_pcsat, id_pcocsat, id_pcq ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
48	!$OMP THREADPRIVATE(id_pcsat, id_pcocsat, id_pcq)
49	INTEGER,SAVE :: id_pcs0, id_pcos0, id_pcq0 ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
50	! ! qui ne sont pas transportes par la convection
51	!$OMP THREADPRIVATE(id_pcs0, id_pcos0, id_pcq0)
52
53	INTEGER, SAVE:: id_o3
54	!$OMP THREADPRIVATE(id_o3)
55	! index of ozone tracer with Cariolle parameterization
56	! 0 means no ozone tracer
57
58	LOGICAL,SAVE :: rnpb=.FALSE. ! Presence du couple Rn222, Pb210
59	!$OMP THREADPRIVATE(rnpb)
60
61
62	CONTAINS
63
64	SUBROUTINE traclmdz_from_restart(trs_in)
65	! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc).
66	! This subroutine is called from phyetat0 after the field trs_in has been read.
67
68	USE dimphy
69	USE infotrac_phy, ONLY: nbtr
70
71	! Input argument
72	REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in
73
74	! Local variables
75	INTEGER :: ierr
76
77	! Allocate restart variables trs
78	ALLOCATE( trs(klon,nbtr), stat=ierr)
79	IF (ierr /= 0) CALL abort_physic('traclmdz_from_restart', 'pb in allocation 1',1)
80
81	! Initialize trs with values read from restart file
82	trs(:,:) = trs_in(:,:)
83
84	END SUBROUTINE traclmdz_from_restart
85
86
87	SUBROUTINE traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
88	! This subroutine allocates and initialize module variables and control variables.
89	! Initialization of the tracers should be done here only for those not found in the restart file.
90	USE dimphy
91	USE infotrac_phy, ONLY: nbtr, nqtot, tracers, pbl_flg, conv_flg
92	USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
93	USE press_coefoz_m, ONLY: press_coefoz
94	USE mod_grid_phy_lmdz
95	USE mod_phys_lmdz_para
96	USE indice_sol_mod
97	USE print_control_mod, ONLY: lunout
98	USE strings_mod, ONLY: strLower
99
100	! Input variables
101	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
102	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes en degres pour chaque point
103	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes en degres pour chaque point
104	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
105	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA]
106	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
107	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
108	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
109	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
110
111	! Output variables
112	LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
113	LOGICAL,INTENT(OUT) :: lessivage
114
115	! Local variables
116	INTEGER :: ierr, it, iq, i, k
117	REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global
118	REAL, DIMENSION(klev) :: mintmp, maxtmp
119	LOGICAL :: zero
120	! RomP >>> profil initial Be7
121	integer ilesfil
122	parameter (ilesfil=1)
123	integer irr,kradio
124	real beryllium(klon,klev)
125	! profil initial Pb210
126	integer ilesfil2
127	parameter (ilesfil2=1)
128	integer irr2,kradio2
129	real plomb(klon,klev)
130	!! RomP <<<
131	! --------------------------------------------
132	! Allocation
133	! --------------------------------------------
134	ALLOCATE( scavtr(nbtr), stat=ierr)
135	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 9',1)
136	scavtr(:)=1.
137
138	ALLOCATE( radio(nbtr), stat=ierr)
139	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 11',1)
140	radio(:) = .false. ! Par defaut pas decroissance radioactive
141
142	ALLOCATE( masktr(klon,nbtr), stat=ierr)
143	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 2',1)
144
145	ALLOCATE( fshtr(klon,nbtr), stat=ierr)
146	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 3',1)
147
148	ALLOCATE( hsoltr(nbtr), stat=ierr)
149	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 4',1)
150
151	ALLOCATE( tautr(nbtr), stat=ierr)
152	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 5',1)
153	tautr(:) = 0.
154
155	ALLOCATE( vdeptr(nbtr), stat=ierr)
156	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 6',1)
157	vdeptr(:) = 0.
158
159
160	lessivage = .TRUE.
161	!!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
162	!! call getin('lessivage',lessivage)
163	!! if(lessivage) then
164	!! print*,'lessivage lsc ON'
165	!! else
166	!! print*,'lessivage lsc OFF'
167	!! endif
168	aerosol(:) = .FALSE. ! Tous les traceurs sont des gaz par defaut
169
170	!
171	! Recherche des traceurs connus : Be7, O3, CO2,...
172	! --------------------------------------------
173	id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
174	id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
175	it = 0
176	DO iq = 1, nqtot
177	IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
178	it = it+1
179	SELECT CASE(strLower(tracers(iq)%name))
180	CASE("rn"); id_rn = it ! radon
181	CASE("pb"); id_pb = it ! plomb
182	CASE("aga"); id_aga = it ! Age of stratospheric air
183	CASE("be","be7");id_be = it ! Recherche du Beryllium 7
184	CASE("o3"); id_o3 = it ! Recherche de l'ozone
185	CASE("pcsat"); id_pcsat = it
186	CASE("pcocsat"); id_pcocsat= it
187	CASE("pcq"); id_pcq = it
188	CASE("pcs0"); id_pcs0 = it
189	CASE("pcos0"); id_pcos0 = it
190	CASE("pcq0"); id_pcq0 = it
191	CASE DEFAULT
192	WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iq)%name)
193	END SELECT
194
195	SELECT CASE(strLower(tracers(iq)%name))
196	CASE("pb") !--- RomP >>> profil initial de PB210
197	OPEN(ilesfil2,file='prof.pb210',status='old',iostat=irr2)
198	IF(irr2 == 0) THEN
199	READ(ilesfil2,*) kradio2
200	WRITE(lunout,*)'number of levels for pb210 profile ',kradio2
201	DO k=kradio2,1,-1; READ (ilesfil2,*) plomb(:,k); END DO
202	CLOSE(ilesfil2)
203	tr_seri(:,:,id_pb) = plomb(:,:)
204	ELSE
205	WRITE(lunout,*)'Prof. Pb210 does not exist: use restart values'
206	END IF
207	CASE("aga")
208	radio(id_aga) = .FALSE.
209	aerosol(id_aga) = .FALSE.
210	pbl_flg(id_aga) = 0
211	! Find the first model layer above 1.5km from the surface
212	IF (klev>=30) THEN
213	lev_1p5km=6 !--- NB: This value is for klev=39
214	ELSE IF (klev>=10) THEN
215	lev_1p5km=5 !--- NB: This value is for klev=19
216	ELSE
217	lev_1p5km=klev/2
218	END IF
219	CASE("be","be7")
220	ALLOCATE( srcbe(klon,klev) )
221	radio(id_be) = .TRUE.
222	aerosol(id_be) = .TRUE. !--- Le Be est un aerosol
223	CALL init_be(pctsrf,pplay,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
224	WRITE(lunout,*) 'Initialisation srcBe: OK'
225	!--- RomP >>> profil initial de Be7
226	OPEN(ilesfil,file='prof.be7',status='old',iostat=irr)
227	IF(irr == 0) THEN
228	READ(ilesfil,*) kradio
229	WRITE(lunout,*)'number of levels for Be7 profile ',kradio
230	DO k=kradio,1,-1; READ(ilesfil,*) beryllium(:,k); END DO
231	CLOSE(ilesfil)
232	tr_seri(:,:,id_be)=beryllium(:,:)
233	ELSE
234	WRITE(lunout,*)'Prof. Be7 does not exist: use restart values'
235	END IF
236	CASE("o3") !--- Parametrisation par la chimie de Cariolle
237	CALL alloc_coefoz !--- Allocate ozone coefficients
238	CALL press_coefoz !--- Read input pressure levels
239	CASE("pcs0","pcos0","pcq0")
240	conv_flg(it)=0 !--- No transport by convection for this tracer
241	END SELECT
242	END DO
243
244	!
245	! Valeurs specifiques pour les traceurs Rn222 et Pb210
246	! ----------------------------------------------
247	IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
248	rnpb = .TRUE.
249	radio(id_rn)= .TRUE.
250	radio(id_pb)= .TRUE.
251	pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
252	pbl_flg(id_pb) = 0 ! au lieu de clsol=true
253	aerosol(id_rn) = .FALSE.
254	aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
255
256	CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
257	END IF
258
259	!
260	! Check if all tracers have restart values
261	! ----------------------------------------------
262	it = 0
263	DO iq = 1, nqtot
264	IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
265	it = it+1
266	! Test if tracer is zero everywhere.
267	! Done by master process MPI and master thread OpenMP
268	CALL gather(tr_seri(:,:,it),varglo)
269	IF (is_mpi_root .AND. is_omp_root) THEN
270	mintmp=MINVAL(varglo,dim=1)
271	maxtmp=MAXVAL(varglo,dim=1)
272	IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
273	! Tracer is zero everywhere
274	zero=.TRUE.
275	ELSE
276	zero=.FALSE.
277	END IF
278	END IF
279
280	! Distribute variable at all processes
281	CALL bcast(zero)
282
283	! Initalize tracer that was not found in restart file.
284	IF (zero) THEN
285	! The tracer was not found in restart file or it was equal zero everywhere.
286	WRITE(lunout,*) "The tracer ",trim(tracers(iq)%name)," will be initialized"
287	IF (it==id_pcsat .OR. it==id_pcq .OR. &
288	it==id_pcs0 .OR. it==id_pcq0) THEN
289	tr_seri(:,:,it) = 100.
290	ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
291	DO i = 1, klon
292	IF ( pctsrf (i, is_oce) == 0. ) THEN
293	tr_seri(i,:,it) = 0.
294	ELSE
295	tr_seri(i,:,it) = 100.
296	END IF
297	END DO
298	ELSE
299	! No specific initialization exist for this tracer
300	tr_seri(:,:,it) = 0.
301	END IF
302	END IF
303	END DO
304
305	END SUBROUTINE traclmdz_init
306
307	SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
308	cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
309	rh, pphi, ustar, wstar, ale_bl, ale_wake, zu10m, zv10m, &
310	tr_seri, source, d_tr_cl,d_tr_dec, zmasse) !RomP
311
312	USE dimphy
313	USE infotrac_phy, ONLY: nbtr, pbl_flg, solsym
314	USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
315	USE o3_chem_m, ONLY: o3_chem
316	USE indice_sol_mod
317
318	INCLUDE "YOMCST.h"
319
320	!==========================================================================
321	! -- DESCRIPTION DES ARGUMENTS --
322	!==========================================================================
323
324	! Input arguments
325	!
326	!Configuration grille,temps:
327	INTEGER,INTENT(IN) :: nstep ! nombre d'appels de la physiq
328	INTEGER,INTENT(IN) :: julien ! Jour julien
329	REAL,INTENT(IN) :: gmtime
330	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
331	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
332	REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
333
334	!
335	!Physique:
336	!--------
337	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
338	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
339	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
340	REAL,intent(in):: zmasse (:, :) ! dim(klon,klev) density of air, in kg/m2
341
342
343	!Couche limite:
344	!--------------
345	!
346	REAL,DIMENSION(klon),INTENT(IN) :: cdragh ! coeff drag pour T et Q
347	REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh ! coeff melange CL (m**2/s)
348	REAL,DIMENSION(klon),INTENT(IN) :: yu1 ! vents au premier niveau
349	REAL,DIMENSION(klon),INTENT(IN) :: yv1 ! vents au premier niveau
350	LOGICAL,INTENT(IN) :: couchelimite
351	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
352	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! Humidite relative
353	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentie
354	REAL,DIMENSION(klon),INTENT(IN) :: ustar ! ustar (m/s)
355	REAL,DIMENSION(klon),INTENT(IN) :: wstar,ale_bl,ale_wake ! wstar (m/s) and Avail. Lifti. Energ.
356	REAL,DIMENSION(klon),INTENT(IN) :: zu10m ! vent zonal 10m (m/s)
357	REAL,DIMENSION(klon),INTENT(IN) :: zv10m ! vent zonal 10m (m/s)
358
359	! Arguments necessaires pour les sources et puits de traceur:
360	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
361	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
362
363	! InOutput argument
364	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
365
366	! Output argument
367	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: source ! a voir lorsque le flux de surface est prescrit
368	REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT) :: d_tr_cl ! Td couche limite/traceur
369
370	!=======================================================================================
371	! -- VARIABLES LOCALES TRACEURS --
372	!=======================================================================================
373
374	INTEGER :: i, k, it
375	INTEGER :: lmt_pas ! number of time steps of "physics" per day
376
377	REAL,DIMENSION(klon) :: d_trs ! Td dans le reservoir
378	REAL,DIMENSION(klon,klev) :: qsat ! pression de la vapeur a saturation
379	REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
380	REAL :: zrho ! Masse Volumique de l'air KgA/m3
381	REAL :: amn, amx
382	!
383	!=================================================================
384	! Ajout de la production en Be7 (Beryllium) srcbe U/s/kgA
385	!=================================================================
386	!
387	IF ( id_be /= 0 ) THEN
388	DO k = 1, klev
389	DO i = 1, klon
390	tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
391	END DO
392	END DO
393	WRITE(,) 'Ajout srcBe dans tr_seri: OK'
394	END IF
395
396
397	!=================================================================
398	! Update pseudo-vapor tracers
399	!=================================================================
400
401	CALL q_sat(klon*klev,t_seri,pplay,qsat)
402
403	IF ( id_pcsat /= 0 ) THEN
404	DO k = 1, klev
405	DO i = 1, klon
406	IF ( pplay(i,k).GE.85000.) THEN
407	tr_seri(i,k,id_pcsat) = qsat(i,k)
408	ELSE
409	tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))
410	END IF
411	END DO
412	END DO
413	END IF
414
415	IF ( id_pcocsat /= 0 ) THEN
416	DO k = 1, klev
417	DO i = 1, klon
418	IF ( pplay(i,k).GE.85000.) THEN
419	IF ( pctsrf (i, is_oce) > 0. ) THEN
420	tr_seri(i,k,id_pcocsat) = qsat(i,k)
421	ELSE
422	tr_seri(i,k,id_pcocsat) = 0.
423	END IF
424	ELSE
425	tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
426	END IF
427	END DO
428	END DO
429	END IF
430
431	IF ( id_pcq /= 0 ) THEN
432	DO k = 1, klev
433	DO i = 1, klon
434	IF ( pplay(i,k).GE.85000.) THEN
435	tr_seri(i,k,id_pcq) = sh(i,k)
436	ELSE
437	tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))
438	END IF
439	END DO
440	END DO
441	END IF
442
443
444	IF ( id_pcs0 /= 0 ) THEN
445	DO k = 1, klev
446	DO i = 1, klon
447	IF ( pplay(i,k).GE.85000.) THEN
448	tr_seri(i,k,id_pcs0) = qsat(i,k)
449	ELSE
450	tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))
451	END IF
452	END DO
453	END DO
454	END IF
455
456
457	IF ( id_pcos0 /= 0 ) THEN
458	DO k = 1, klev
459	DO i = 1, klon
460	IF ( pplay(i,k).GE.85000.) THEN
461	IF ( pctsrf (i, is_oce) > 0. ) THEN
462	tr_seri(i,k,id_pcos0) = qsat(i,k)
463	ELSE
464	tr_seri(i,k,id_pcos0) = 0.
465	END IF
466	ELSE
467	tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
468	END IF
469	END DO
470	END DO
471	END IF
472
473
474	IF ( id_pcq0 /= 0 ) THEN
475	DO k = 1, klev
476	DO i = 1, klon
477	IF ( pplay(i,k).GE.85000.) THEN
478	tr_seri(i,k,id_pcq0) = sh(i,k)
479	ELSE
480	tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
481	END IF
482	END DO
483	END DO
484	END IF
485
486	!=================================================================
487	! Update tracer : Age of stratospheric air
488	!=================================================================
489	IF (id_aga/=0) THEN
490
491	! Bottom layers
492	DO k = 1, lev_1p5km
493	tr_seri(:,k,id_aga) = 0.0
494	END DO
495
496	! Layers above 1.5km
497	DO k = lev_1p5km+1,klev-1
498	tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
499	END DO
500
501	! Top layer
502	tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
503
504	END IF
505
506	!======================================================================
507	! -- Calcul de l'effet de la couche limite --
508	!======================================================================
509
510	IF (couchelimite) THEN
511	source(:,:) = 0.0
512
513	IF (id_be /=0) THEN
514	DO i=1, klon
515	zrho = pplay(i,1)/t_seri(i,1)/RD
516	source(i,id_be) = - vdeptr(id_be)tr_seri(i,1,id_be)zrho
517	END DO
518	END IF
519
520	END IF
521
522	DO it=1, nbtr
523	IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
524	! couche limite avec quantite dans le sol calculee
525	CALL cltracrn(it, pdtphys, yu1, yv1, &
526	cdragh, coefh,t_seri,ftsol,pctsrf, &
527	tr_seri(:,:,it),trs(:,it), &
528	paprs, pplay, zmasse * rg, &
529	masktr(:,it),fshtr(:,it),hsoltr(it),&
530	tautr(it),vdeptr(it), &
531	xlat,d_tr_cl(:,:,it),d_trs)
532
533	DO k = 1, klev
534	DO i = 1, klon
535	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
536	END DO
537	END DO
538
539	! Traceur dans le reservoir sol
540	DO i = 1, klon
541	trs(i,it) = trs(i,it) + d_trs(i)
542	END DO
543	END IF
544	END DO
545
546
547	!======================================================================
548	! Calcul de l'effet du puits radioactif
549	!======================================================================
550	CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
551
552	DO it=1,nbtr
553	WRITE(solsym(it),'(i2)') it
554	END DO
555
556	DO it=1,nbtr
557	IF(radio(it)) then
558	DO k = 1, klev
559	DO i = 1, klon
560	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
561	END DO
562	END DO
563	CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
564	END IF
565	END DO
566
567	!======================================================================
568	! Parameterization of ozone chemistry
569	!======================================================================
570
571	IF (id_o3 /= 0) then
572	lmt_pas = NINT(86400./pdtphys)
573	IF (MOD(nstep - 1, lmt_pas) == 0) THEN
574	! Once per day, update the coefficients for ozone chemistry:
575	CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
576	END IF
577	CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
578	xlon, tr_seri(:, :, id_o3))
579	END IF
580
581	END SUBROUTINE traclmdz
582
583
584	SUBROUTINE traclmdz_to_restart(trs_out)
585	! This subroutine is called from phyredem.F where the module
586	! variable trs is written to restart file (restartphy.nc)
587	USE dimphy
588	USE infotrac_phy, ONLY: nbtr
589
590	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
591	INTEGER :: ierr
592
593	IF ( ALLOCATED(trs) ) THEN
594	trs_out(:,:) = trs(:,:)
595	ELSE
596	! No previous allocate of trs. This is the case for create_etat0_limit.
597	trs_out(:,:) = 0.0
598	END IF
599
600	END SUBROUTINE traclmdz_to_restart
601
602
603	END MODULE traclmdz_mod

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

Download in other formats:

Original Format