Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

set_ub_vals.F90 @ 6610

Last change on this file since 6610 was 6610, checked in by acosce, 10 months ago
INCA used for ICOLMDZORINCA_CO2_Transport_GMD_2023
File size: 61.1 KB

Line
1	!$Id: set_ub_vals.F90 157 2010-01-18 14:21:59Z acosce $
2	!! =========================================================================
3	!! INCA - INteraction with Chemistry and Aerosols
4	!!
5	!! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE)
6	!! Unite mixte CEA-CNRS-UVSQ
7	!!
8	!! Contributors to this INCA subroutine:
9	!!
10	!! Didier Hauglustaine, LSCE, hauglustaine@cea.fr
11	!!
12	!! Anne Cozic, LSCE, anne.cozic@cea.fr
13	!! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr
14	!!
15	!! This software is a computer program whose purpose is to simulate the
16	!! atmospheric gas phase and aerosol composition. The model is designed to be
17	!! used within a transport model or a general circulation model. This version
18	!! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts
19	!! for emissions, transport (resolved and sub-grid scale), photochemical
20	!! transformations, and scavenging (dry deposition and washout) of chemical
21	!! species and aerosols interactively in the GCM. Several versions of the INCA
22	!! model are currently used depending on the envisaged applications with the
23	!! chemistry-climate model.
24	!!
25	!! This software is governed by the CeCILL license under French law and
26	!! abiding by the rules of distribution of free software. You can use,
27	!! modify and/ or redistribute the software under the terms of the CeCILL
28	!! license as circulated by CEA, CNRS and INRIA at the following URL
29	!! "http://www.cecill.info".
30	!!
31	!! As a counterpart to the access to the source code and rights to copy,
32	!! modify and redistribute granted by the license, users are provided only
33	!! with a limited warranty and the software's author, the holder of the
34	!! economic rights, and the successive licensors have only limited
35	!! liability.
36	!!
37	!! In this respect, the user's attention is drawn to the risks associated
38	!! with loading, using, modifying and/or developing or reproducing the
39	!! software by the user in light of its specific status of free software,
40	!! that may mean that it is complicated to manipulate, and that also
41	!! therefore means that it is reserved for developers and experienced
42	!! professionals having in-depth computer knowledge. Users are therefore
43	!! encouraged to load and test the software's suitability as regards their
44	!! requirements in conditions enabling the security of their systems and/or
45	!! data to be ensured and, more generally, to use and operate it in the
46	!! same conditions as regards security.
47	!!
48	!! The fact that you are presently reading this means that you have had
49	!! knowledge of the CeCILL license and that you accept its terms.
50	!! =========================================================================
51
52	#include <inca_define.h>
53
54	#if defined(AERONLY) \|\| defined(GES)
55
56
57
58
59	SUBROUTINE XIOS_OXYDANT_READ(calday)
60	!----------------------------------------------------------------------
61	! ... Read OXYDANT distributions
62	! Didier Hauglustaine, IPSL, 1999.
63	!----------------------------------------------------------------------
64
65	USE OXYDANT_COM
66	USE MOD_GRID_INCA
67	USE INCA_DIM
68	USE MOD_INCA_PARA
69	USE PRINT_INCA
70	USE XIOS_INCA
71	IMPLICIT NONE
72
73	INCLUDE 'netcdf.inc'
74
75	!----------------------------------------------------------------------
76	! ... Dummy args
77	!----------------------------------------------------------------------
78	REAL, INTENT(IN) :: calday
79
80	!----------------------------------------------------------------------
81	! ... Local variables
82	!----------------------------------------------------------------------
83	INTEGER :: iret
84	INTEGER :: varid
85	INTEGER :: oldlotime, oldhitime
86	INTEGER, DIMENSION(4) :: start3, count3
87	LOGICAL, SAVE :: first = .TRUE.
88	!$OMP THREADPRIVATE(first)
89
90	REAL, DIMENSION(12) :: timeo_inter
91	REAL :: dt, dt1, tint
92	INTEGER :: lonlev
93
94
95	IF (first) THEN
96
97
98
99	ALLOCATE(timeo(ntime_oxyd))
100	ALLOCATE(o3oxyd_year(PLON, PLEV, ntime_oxyd))
101	ALLOCATE(o1doxyd_year(PLON, PLEV, ntime_oxyd))
102	ALLOCATE(hno3oxyd_year(PLON, PLEV, ntime_oxyd))
103	ALLOCATE(ohoxyd_year(PLON, PLEV, ntime_oxyd))
104	ALLOCATE(h2o2oxyd_year(PLON, PLEV, ntime_oxyd))
105	ALLOCATE(no2oxyd_year(PLON, PLEV, ntime_oxyd))
106	ALLOCATE(no3oxyd_year(PLON, PLEV, ntime_oxyd))
107
108	! write(lunout,*) '--------------'
109	! write(lunout,*) 'lecture oxyd : '
110	! write(lunout,*) 'PLON = ', PLON
111	! write(lunout,*) 'PLEV = ', PLEV
112	! write(lunout,*) 'ntime_oxy = ', ntime_oxyd
113	! write(lunout,*) '--------------'
114
115	CALL xios_inca_change_context("inca")
116
117	CALL xios_inca_recv_field_glo("oxydtime_read", timeo )
118
119	CALL xios_inca_recv_field("O3_interp" , o3oxyd_year )
120	CALL xios_inca_recv_field("O1D_interp" , o1doxyd_year )
121	CALL xios_inca_recv_field("HNO3_interp" , hno3oxyd_year)
122	CALL xios_inca_recv_field("OH_interp" , ohoxyd_year )
123	CALL xios_inca_recv_field("H2O2_interp" , h2o2oxyd_year)
124	CALL xios_inca_recv_field("NO2_interp" , no2oxyd_year )
125	CALL xios_inca_recv_field("NO3_interp" , no3oxyd_year )
126
127	CALL xios_inca_send_field("NO3_oxyd_read", no3oxyd_year(:,:,12))
128
129	CALL xios_inca_change_context("LMDZ")
130
131
132
133	first = .false.
134	ENDIF
135
136
137	lonlev = PLON * nlevo
138
139
140	! ... Check to see if model time still bounded by data set
141	oldlotime = lotime
142	oldhitime = hitime
143	timeo_inter(:) = mod(timeo(:)/86400, 365.)
144
145
146	CALL FINDPLB (timeo_inter, ntime_oxyd, calday, lotime)
147
148	IF ( cyclical ) THEN
149	hitime = MOD(lotime,ntime_oxyd) + 1
150	ELSE
151	hitime = lotime + 1
152	END IF
153	! -------------------- interpolation sur le temps
154	!----------------------------------------------------------------------
155	! Note: 365 day year inconsistent with LMDz (360 days) !!!
156	!----------------------------------------------------------------------
157	! ... First interpolate linearly on current model time
158	!
159	!
160	IF ( timeo(hitime) < timeo(lotime) ) THEN
161	dt = 365. - mod(timeo(lotime)/(86400),365.) + mod(timeo(hitime)/(86400),365.)
162	IF (calday <= mod(timeo(hitime)/(86400),365.) ) THEN
163	dt1 = 365. - mod(timeo(lotime)/(86400),365.) + calday
164	ELSE
165	dt1 = calday - mod(timeo(lotime)/(86400),365.)
166	END IF
167	ELSE
168	dt = mod(timeo(hitime)/(86400),365.) - mod(timeo(lotime)/(86400),365.)
169	dt1 = calday - mod(timeo(lotime)/(86400),365.)
170	END IF
171
172	tint = dt1/dt
173
174	ohoxyd(:,:) = ohoxyd_year(:,:,lotime) + (ohoxyd_year(:,:,hitime)-ohoxyd_year(:,:,lotime))*tint
175	o1doxyd(:,:) = o1doxyd_year(:,:,lotime) + (o1doxyd_year(:,:,hitime)-o1doxyd_year(:,:,lotime))*tint
176	o3oxyd(:,:) = o3oxyd_year(:,:,lotime) + (o3oxyd_year(:,:,hitime)-o3oxyd_year(:,:,lotime))*tint
177	no3oxyd(:,:) = no3oxyd_year(:,:,lotime) + (no3oxyd_year(:,:,hitime)-no3oxyd_year(:,:,lotime))*tint
178	h2o2oxyd(:,:) = h2o2oxyd_year(:,:,lotime) + (h2o2oxyd_year(:,:,hitime)-h2o2oxyd_year(:,:,lotime))*tint
179	hno3oxyd(:,:) = hno3oxyd_year(:,:,lotime) + (hno3oxyd_year(:,:,hitime)-hno3oxyd_year(:,:,lotime))*tint
180	no2oxyd(:,:) = no2oxyd_year(:,:,lotime) + (no2oxyd_year(:,:,hitime)-no2oxyd_year(:,:,lotime))*tint
181
182	CALL xios_inca_change_context("inca")
183	CALL xios_inca_send_field("O3_oxyd" , o3oxyd )
184	CALL xios_inca_send_field("O1D_oxyd" , o1doxyd )
185	CALL xios_inca_send_field("HNO3_oxyd" , hno3oxyd )
186	CALL xios_inca_send_field("OH_oxyd" , ohoxyd )
187	CALL xios_inca_send_field("H2O2_oxyd" , h2o2oxyd )
188	CALL xios_inca_send_field("NO2_oxyd" , no2oxyd )
189	CALL xios_inca_send_field("NO3_oxyd" , no3oxyd )
190	CALL xios_inca_change_context("LMDZ")
191
192
193	END SUBROUTINE XIOS_OXYDANT_READ
194
195
196	#else
197	! IF not defined aeronly and ges
198
199	SUBROUTINE OZCLIM_READ(calday)
200	!----------------------------------------------------------------------
201	! ... Read ozone distributions
202	! Didier Hauglustaine, IPSL, 1999.
203	!----------------------------------------------------------------------
204
205	USE O3CLIM_COM
206	USE INCA_DIM
207	USE MOD_INCA_PARA
208	USE PRINT_INCA
209
210	IMPLICIT NONE
211
212	INCLUDE 'netcdf.inc'
213
214	!----------------------------------------------------------------------
215	! ... Dummy args
216	!----------------------------------------------------------------------
217	REAL, INTENT(IN) :: calday
218
219	!----------------------------------------------------------------------
220	! ... Local variables
221	!----------------------------------------------------------------------
222	INTEGER :: iret
223	INTEGER :: varid
224	INTEGER :: oldlotime, oldhitime
225	INTEGER, DIMENSION(3) :: start, count
226	REAL :: o3climbd_glo(nbp_glo,nlevc,2)
227
228	! ... Check to see if model time still bounded by data set
229	oldlotime = lotime
230	oldhitime = hitime
231	CALL FINDPLB (timec, ntimec, calday, lotime)
232	IF ( cyclical ) THEN
233	hitime = MOD(lotime,ntimec) + 1
234	ELSE
235	hitime = lotime + 1
236	END IF
237
238	! ... Read new data if needed
239	IF ( hitime /= oldhitime ) THEN
240	loind = hiind
241	hiind = MOD( loind, 2) + 1
242
243	!$OMP MASTER
244	IF (is_mpi_root) THEN
245	iret = nf_inq_varid(ncidc, 'O3', varid)
246	CALL check_err(iret, 'OZCLIM_READ','problem to check varid O3')
247
248	count(1) = nbp_glo
249	count(2) = nlevc
250	count(3) = 1
251	start(1) = 1
252	start(2) = 1
253
254	WRITE(lunout,*) 'OZCLIM_READ : read new data for time ', timec(hitime)
255	start(3) = hitime
256	iret = nf_get_vara_double(ncidc, varid, start, count,o3climbd_glo(1,1,hiind))
257	CALL check_err(iret, 'OZCLIM_READ', 'problem to read O3 hiind')
258	ENDIF
259	!$OMP END MASTER
260	CALL scatter(o3climbd_glo(:,:,hiind),o3climbd(:,:,hiind))
261
262
263	IF ( lotime /= oldhitime ) THEN
264	!$OMP MASTER
265	IF (is_mpi_root) THEN
266	WRITE(lunout,*) 'OZCLIM_READ : read new data for time ', timec(lotime)
267	start(3) = lotime
268	iret = nf_get_vara_double(ncidc, varid, start, count, o3climbd_glo(1,1,loind))
269	CALL check_err(iret, 'OZCLIM_READ', 'problem to read O3 loind')
270
271	ENDIF
272	!$OMP END MASTER
273	CALL scatter(o3climbd_glo(:,:,loind),o3climbd(:,:,loind))
274
275	END IF
276
277	END IF
278
279	END SUBROUTINE OZCLIM_READ
280
281	SUBROUTINE OZCLIM_INTERP(calday, pmid)
282	!----------------------------------------------------------------------
283	! ... Interpolate ozone on model time and on pressure levels
284	! Didier Hauglustaine, IPSL, 1999.
285	!----------------------------------------------------------------------
286
287	USE O3CLIM_COM
288	USE INCA_DIM
289	IMPLICIT NONE
290
291	!----------------------------------------------------------------------
292	! ... Dummy args
293	!----------------------------------------------------------------------
294	REAL, INTENT(IN) :: calday
295	REAL, INTENT(IN) :: pmid(PLON,PLEV)
296
297	!----------------------------------------------------------------------
298	! ... Local variables
299	!----------------------------------------------------------------------
300	REAL :: dt, dt1, tint
301	REAL :: rma = 48./28.
302	REAL :: ploc(klonc,nlevc)
303	INTEGER :: lonlev
304	INTEGER :: i
305	INTEGER :: index(PLON,PLEV)
306
307	index = -9999
308	lonlev = klonc * nlevc
309
310	DO i = 1, klonc
311	ploc(i,:) = presc(:)
312	END DO
313
314	!----------------------------------------------------------------------
315	! Note: 365 day year inconsistent with LMDz (360 days) !!!
316	!----------------------------------------------------------------------
317	! ... First interpolate linearly on current model time
318
319	IF ( timec(hitime) < timec(lotime) ) THEN
320	dt = 365. - timec(lotime) + timec(hitime)
321	IF (calday <= timec(hitime) ) THEN
322	dt1 = 365. - timec(lotime) + calday
323	ELSE
324	dt1 = calday - timec(lotime)
325	END IF
326	ELSE
327	dt = timec(hitime) - timec(lotime)
328	dt1 = calday - timec(lotime)
329	END IF
330
331	tint = dt1/dt
332
333	CALL linintp (lonlev, &
334	0., &
335	1., &
336	tint, &
337	o3climbd(1,1,loind), &
338	o3climbd(1,1,hiind), &
339	o3climcr)
340
341	! ... Put on model pressure levels
342
343	CALL pinterp (PLON, &
344	o3climcr, &
345	ploc, &
346	nlevc, &
347	o3clim, &
348	pmid, &
349	PLEV, &
350	index, &
351	.false. )
352
353	! ... vmr to mmr
354
355	o3clim = o3clim * rma
356
357	END SUBROUTINE OZCLIM_INTERP
358
359
360	SUBROUTINE OZRESET (mmr, pmid, temp)
361	!----------------------------------------------------------------------
362	! ... Reset O3 boundary conditions for use in radiation
363	! Didier Hauglustaine, IPSL, 2002.
364	!----------------------------------------------------------------------
365
366	USE O3CLIM_COM
367	USE SPECIES_NAMES
368	USE OXYDANT_COM
369	USE AIRPLANE_SRC, ONLY : itrop
370	USE CHEM_MODS, ONLY : adv_mass
371	USE INCA_DIM
372	IMPLICIT NONE
373
374	!----------------------------------------------------------------------
375	! ... Dummy args
376	!----------------------------------------------------------------------
377	REAL, INTENT(IN) :: pmid(PLON,PLEV)
378	REAL, INTENT(IN) :: temp(PLON,PLEV)
379	REAL, INTENT(INOUT) :: mmr(PLON,PLEV,PCNST)
380
381	!----------------------------------------------------------------------
382	! ... Local variables
383	!----------------------------------------------------------------------
384	INTEGER :: i, k
385	INTEGER :: ktrop, krelax, k380
386	REAL, PARAMETER :: dry_mass = 28.966
387	REAL :: theta(PLON,PLEV)
388
389	theta(:,:) = temp(:,:)(1.e5/pmid(:,:))*0.286
390
391	DO i = 1, PLON
392
393	ktrop = MAX(1,itrop(i))
394
395	DO k = ktrop, PLEV
396	IF (theta(i,k) >= 380.) THEN
397	k380 = k
398	EXIT
399	ENDIF
400	ENDDO
401
402	! krelax = MIN(ktrop+2,PLEV) ! 2 levels above tropopause
403	krelax = k380 ! 380K theta surface
404
405	mmr(i,krelax:PLEV,id_O3) = o3clim(i,krelax:PLEV)
406	END DO
407
408	END SUBROUTINE OZRESET
409
410	SUBROUTINE FIELD_PREP (mmr)
411	!----------------------------------------------------------------------
412	! ... Prepare fields to be used in LMDz radiation
413	! Didier Hauglustaine, IPSL, 2002. 2018.
414	!----------------------------------------------------------------------
415
416	USE SPECIES_NAMES
417	USE CHEM_MODS, ONLY : o3_inca, h2o_inca, ch4_inca, n2o_inca, cfc11_inca, cfc12_inca, &
418	ch4_inca_2d, n2o_inca_2d, cfc11_inca_2d, cfc12_inca_2d, adv_mass
419	USE INCA_DIM
420	IMPLICIT NONE
421
422	!----------------------------------------------------------------------
423	! ... Dummy args
424	!----------------------------------------------------------------------
425	REAL, INTENT(IN) :: mmr(PLON,PLEV,PCNST)
426
427	!----------------------------------------------------------------------
428	! ... Local variables
429	!----------------------------------------------------------------------
430	INTEGER :: jj, jk, jki, jkf
431	INTEGER, PARAMETER :: ng1 = 2, ng1p1 = ng1 + 1 !From raddimlw.h
432
433	! Fields to be used in LMDz in mass mixing ratio
434	#ifdef NMHC
435	o3_inca(:,:) = mmr(:,:,id_O3)
436	#endif
437	! ch4_inca(:,:) = mmr(:,:,id_CH4)
438	! n2o_inca(:,:) = mmr(:,:,id_N2O)
439	! cfc11_inca(:,:) = mmr(:,:,id_CFC11)
440	! cfc12_inca(:,:) = mmr(:,:,id_CFC12)
441	! h2o_inca(:,:) = mmr(:,:,id_H2O)
442	!
443	!
444	!
445	! ! Ozone field POZON (kg/kg)
446	! ! Note: For POZON (used in radlwsw) simply use o3_inca from above (in mmr)
447	!
448	! ! 2D fields for CH4, N2O, CFC11 and CFC12 (in mmr)
449	! DO jj = 1 , PLEV
450	! jki=(jj-1)*ng1p1+1
451	! jkf=jki+ng1
452	! DO jk = jki, jkf
453	! ch4_inca_2d(:,jk) = ch4_inca(:,jj)
454	! n2o_inca_2d(:,jk) = n2o_inca(:,jj)
455	! cfc11_inca_2d(:,jk) = cfc11_inca(:,jj)
456	! cfc12_inca_2d(:,jk) = cfc12_inca(:,jj)
457	! ENDDO
458	! ENDDO
459
460	END SUBROUTINE FIELD_PREP
461
462
463	SUBROUTINE OZLIN_READ(calday)
464	!----------------------------------------------------------------------
465	! ... Read ozone linear coefficients
466	! Didier Hauglustaine, IPSL, 1999.
467	!----------------------------------------------------------------------
468
469	USE O3LIN_COM
470	USE INCA_DIM
471	USE MOD_INCA_PARA
472	USE PRINT_INCA
473
474	IMPLICIT NONE
475
476	INCLUDE 'netcdf.inc'
477
478	!----------------------------------------------------------------------
479	! ... Dummy args
480	!----------------------------------------------------------------------
481	REAL, INTENT(IN) :: calday
482
483	!----------------------------------------------------------------------
484	! ... Local variables
485	!----------------------------------------------------------------------
486	INTEGER :: iret
487	INTEGER :: varid1,varid2,varid3,varid4,varid5
488	INTEGER :: varid6,varid7,varid8,varid9
489	INTEGER :: oldlotime, oldhitime
490	INTEGER, DIMENSION(3) :: start, count
491	REAL :: A1bd_glo(nbp_glo,nlevl,2)
492	REAL :: A2bd_glo(nbp_glo,nlevl,2)
493	REAL :: A3bd_glo(nbp_glo,nlevl,2)
494	REAL :: A4bd_glo(nbp_glo,nlevl,2)
495	REAL :: A5bd_glo(nbp_glo,nlevl,2)
496	REAL :: A6bd_glo(nbp_glo,nlevl,2)
497	REAL :: A7bd_glo(nbp_glo,nlevl,2)
498	REAL :: A8bd_glo(nbp_glo,nlevl,2)
499	REAL :: A9bd_glo(nbp_glo,nlevl,2)
500
501	! ... Check to see if model time still bounded by data set
502	oldlotime = lotime
503	oldhitime = hitime
504	CALL FINDPLB (timel, ntimel, calday, lotime)
505	IF ( cyclical ) THEN
506	hitime = MOD(lotime,ntimel) + 1
507	ELSE
508	hitime = lotime + 1
509	END IF
510
511	! ... Read new data if needed
512	IF ( hitime /= oldhitime ) THEN
513	loind = hiind
514	hiind = MOD( loind, 2) + 1
515
516	!$OMP MASTER
517	IF (is_mpi_root) THEN
518	iret = nf_inq_varid(ncidl, 'A1', varid1)
519	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A1')
520	iret = nf_inq_varid(ncidl, 'A2', varid2)
521	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A2')
522	iret = nf_inq_varid(ncidl, 'A3', varid3)
523	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A3')
524	iret = nf_inq_varid(ncidl, 'A4', varid4)
525	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A4')
526	iret = nf_inq_varid(ncidl, 'A5', varid5)
527	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A5')
528	iret = nf_inq_varid(ncidl, 'A6', varid6)
529	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A6')
530	iret = nf_inq_varid(ncidl, 'A7', varid7)
531	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A7')
532	iret = nf_inq_varid(ncidl, 'A8', varid8)
533	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A8')
534	iret = nf_inq_varid(ncidl, 'A9', varid9)
535	CALL check_err(iret, 'OZLIN_READ', 'problem to get id for A9')
536
537	count(1) = nbp_glo
538	count(2) = nlevl
539	count(3) = 1
540	start(1) = 1
541	start(2) = 1
542
543	WRITE(lunout,*) 'OZLIN_READ : read new data for time ', timel(hitime)
544	start(3) = hitime
545	iret = nf_get_vara_double(ncidl, varid1, start, count, A1bd_glo(1,1,hiind))
546	CALL check_err(iret, 'OCLIN_READ', 'problem to read A1bd for hiind')
547	iret = nf_get_vara_double(ncidl, varid2, start, count, A2bd_glo(1,1,hiind))
548	CALL check_err(iret, 'OCLIN_READ', 'problem to read A2bd for hiind')
549	iret = nf_get_vara_double(ncidl, varid3, start, count, A3bd_glo(1,1,hiind))
550	CALL check_err(iret, 'OCLIN_READ', 'problem to read A3bd for hiind')
551	iret = nf_get_vara_double(ncidl, varid4, start, count, A4bd_glo(1,1,hiind))
552	CALL check_err(iret, 'OCLIN_READ', 'problem to read A4bd for hiind')
553	iret = nf_get_vara_double(ncidl, varid5, start, count, A5bd_glo(1,1,hiind))
554	CALL check_err(iret, 'OCLIN_READ', 'problem to read A5bd for hiind')
555	iret = nf_get_vara_double(ncidl, varid6, start, count, A6bd_glo(1,1,hiind))
556	CALL check_err(iret, 'OCLIN_READ', 'problem to read A6bd for hiind')
557	iret = nf_get_vara_double(ncidl, varid7, start, count, A7bd_glo(1,1,hiind))
558	CALL check_err(iret, 'OCLIN_READ', 'problem to read A7bd for hiind')
559	iret = nf_get_vara_double(ncidl, varid8, start, count, A8bd_glo(1,1,hiind))
560	CALL check_err(iret, 'OCLIN_READ', 'problem to read A8bd for hiind')
561	iret = nf_get_vara_double(ncidl, varid9, start, count, A9bd_glo(1,1,hiind))
562	CALL check_err(iret, 'OCLIN_READ', 'problem to read A9bd for hiind')
563
564	ENDIF
565	!$OMP END MASTER
566
567	CALL scatter(A1bd_glo(:,:,hiind),A1bd(:,:,hiind))
568	CALL scatter(A2bd_glo(:,:,hiind),A2bd(:,:,hiind))
569	CALL scatter(A3bd_glo(:,:,hiind),A3bd(:,:,hiind))
570	CALL scatter(A4bd_glo(:,:,hiind),A4bd(:,:,hiind))
571	CALL scatter(A5bd_glo(:,:,hiind),A5bd(:,:,hiind))
572	CALL scatter(A6bd_glo(:,:,hiind),A6bd(:,:,hiind))
573	CALL scatter(A7bd_glo(:,:,hiind),A7bd(:,:,hiind))
574	CALL scatter(A8bd_glo(:,:,hiind),A8bd(:,:,hiind))
575	CALL scatter(A9bd_glo(:,:,hiind),A9bd(:,:,hiind))
576
577	IF ( lotime /= oldhitime ) THEN
578	!$OMP MASTER
579	IF (is_mpi_root) THEN
580	WRITE(lunout,*) 'OZLIN_READ : read new data for time ', timel(lotime)
581	start(3) = lotime
582	iret = nf_get_vara_double(ncidl, varid1, start, count, A1bd_glo(1,1,loind))
583	CALL check_err(iret, 'OZLIN_READ', 'problem to read A1bd for loind')
584	iret = nf_get_vara_double(ncidl, varid2, start, count, A2bd_glo(1,1,loind))
585	CALL check_err(iret, 'OZLIN_READ', 'problem to read A2bd for loind')
586	iret = nf_get_vara_double(ncidl, varid3, start, count, A3bd_glo(1,1,loind))
587	CALL check_err(iret, 'OZLIN_READ', 'problem to read A3bd for loind')
588	iret = nf_get_vara_double(ncidl, varid4, start, count, A4bd_glo(1,1,loind))
589	CALL check_err(iret, 'OZLIN_READ', 'problem to read A4bd for loind')
590	iret = nf_get_vara_double(ncidl, varid5, start, count, A5bd_glo(1,1,loind))
591	CALL check_err(iret, 'OZLIN_READ', 'problem to read A5bd for loind')
592	iret = nf_get_vara_double(ncidl, varid6, start, count, A6bd_glo(1,1,loind))
593	CALL check_err(iret, 'OZLIN_READ', 'problem to read A6bd for loind')
594	iret = nf_get_vara_double(ncidl, varid7, start, count, A7bd_glo(1,1,loind))
595	CALL check_err(iret, 'OZLIN_READ', 'problem to read A7bd for loind')
596	iret = nf_get_vara_double(ncidl, varid8, start, count, A8bd_glo(1,1,loind))
597	CALL check_err(iret, 'OZLIN_READ', 'problem to read A8bd for loind')
598	iret = nf_get_vara_double(ncidl, varid9, start, count, A9bd_glo(1,1,loind))
599	CALL check_err(iret, 'OZLIN_READ', 'problem to read A9bd for loind')
600
601	ENDIF
602	!$OMP END MASTER
603	CALL scatter(A1bd_glo(:,:,loind),A1bd(:,:,loind))
604	CALL scatter(A2bd_glo(:,:,loind),A2bd(:,:,loind))
605	CALL scatter(A3bd_glo(:,:,loind),A3bd(:,:,loind))
606	CALL scatter(A4bd_glo(:,:,loind),A4bd(:,:,loind))
607	CALL scatter(A5bd_glo(:,:,loind),A5bd(:,:,loind))
608	CALL scatter(A6bd_glo(:,:,loind),A6bd(:,:,loind))
609	CALL scatter(A7bd_glo(:,:,loind),A7bd(:,:,loind))
610	CALL scatter(A8bd_glo(:,:,loind),A8bd(:,:,loind))
611	CALL scatter(A9bd_glo(:,:,loind),A9bd(:,:,loind))
612
613	END IF
614
615	END IF
616
617	END SUBROUTINE OZLIN_READ
618
619	SUBROUTINE OZLIN_INTERP(calday, pmid)
620	!----------------------------------------------------------------------
621	! ... Interpolate ozone linear coefficients on model time and on pressure levels
622	!R. Valorso
623	!----------------------------------------------------------------------
624
625	USE O3LIN_COM
626	USE INCA_DIM
627	IMPLICIT NONE
628
629	!----------------------------------------------------------------------
630	! ... Dummy args
631	!----------------------------------------------------------------------
632	REAL, INTENT(IN) :: calday
633	REAL, INTENT(IN) :: pmid(PLON,PLEV)
634
635	!----------------------------------------------------------------------
636	! ... Local variables
637	!----------------------------------------------------------------------
638	REAL :: dt, dt1, tint
639	REAL :: ploc(klonl,nlevl)
640	INTEGER :: lonlev
641	INTEGER :: i
642	INTEGER :: index(PLON,PLEV)
643
644	index = -9999
645	lonlev = klonl * nlevl
646
647	DO i = 1, klonl
648	ploc(i,:) = presl(:)
649	END DO
650
651	!----------------------------------------------------------------------
652	! Note: 365 day year inconsistent with LMDz (360 days) !!!
653	!----------------------------------------------------------------------
654	! ... First interpolate linearly on current model time
655
656	IF ( timel(hitime) < timel(lotime) ) THEN
657	dt = 365. - timel(lotime) + timel(hitime)
658	IF (calday <= timel(hitime) ) THEN
659	dt1 = 365. - timel(lotime) + calday
660	ELSE
661	dt1 = calday - timel(lotime)
662	END IF
663	ELSE
664	dt = timel(hitime) - timel(lotime)
665	dt1 = calday - timel(lotime)
666	END IF
667
668	tint = dt1/dt
669
670	CALL linintp (lonlev, &
671	0., &
672	1., &
673	tint, &
674	A1bd(1,1,loind), &
675	A1bd(1,1,hiind), &
676	A1cr)
677
678	! ... Put on model pressure levels
679
680	CALL pinterp (PLON, &
681	A1cr, &
682	ploc, &
683	nlevl, &
684	A1, &
685	pmid, &
686	PLEV, &
687	index, &
688	.false. )
689
690	! ...
691
692	CALL linintp (lonlev, &
693	0., &
694	1., &
695	tint, &
696	A2bd(1,1,loind), &
697	A2bd(1,1,hiind), &
698	A2cr)
699
700	! ... Put on model pressure levels
701
702	index = -9999
703	CALL pinterp (PLON, &
704	A2cr, &
705	ploc, &
706	nlevl, &
707	A2, &
708	pmid, &
709	PLEV, &
710	index, &
711	.false. )
712	! ...
713
714	CALL linintp (lonlev, &
715	0., &
716	1., &
717	tint, &
718	A3bd(1,1,loind), &
719	A3bd(1,1,hiind), &
720	A3cr)
721
722	! ... Put on model pressure levels
723
724	index = -9999
725	CALL pinterp (PLON, &
726	A3cr, &
727	ploc, &
728	nlevl, &
729	A3, &
730	pmid, &
731	PLEV, &
732	index, &
733	.false. )
734	! ...
735
736	CALL linintp (lonlev, &
737	0., &
738	1., &
739	tint, &
740	A4bd(1,1,loind), &
741	A4bd(1,1,hiind), &
742	A4cr)
743
744	! ... Put on model pressure levels
745
746	index = -9999
747	CALL pinterp (PLON, &
748	A4cr, &
749	ploc, &
750	nlevl, &
751	A4, &
752	pmid, &
753	PLEV, &
754	index, &
755	.false. )
756	! ...
757
758	CALL linintp (lonlev, &
759	0., &
760	1., &
761	tint, &
762	A5bd(1,1,loind), &
763	A5bd(1,1,hiind), &
764	A5cr)
765
766	! ... Put on model pressure levels
767
768	index = -9999
769	CALL pinterp (PLON, &
770	A5cr, &
771	ploc, &
772	nlevl, &
773	A5, &
774	pmid, &
775	PLEV, &
776	index, &
777	.false. )
778	! ...
779
780	CALL linintp (lonlev, &
781	0., &
782	1., &
783	tint, &
784	A6bd(1,1,loind), &
785	A6bd(1,1,hiind), &
786	A6cr)
787
788	! ... Put on model pressure levels
789
790	index = -9999
791	CALL pinterp (PLON, &
792	A6cr, &
793	ploc, &
794	nlevl, &
795	A6, &
796	pmid, &
797	PLEV, &
798	index, &
799	.false. )
800	! ...
801
802	CALL linintp (lonlev, &
803	0., &
804	1., &
805	tint, &
806	A7bd(1,1,loind), &
807	A7bd(1,1,hiind), &
808	A7cr)
809
810	! ... Put on model pressure levels
811
812	index = -9999
813	CALL pinterp (PLON, &
814	A7cr, &
815	ploc, &
816	nlevl, &
817	A7, &
818	pmid, &
819	PLEV, &
820	index, &
821	.false. )
822	! ...
823
824	CALL linintp (lonlev, &
825	0., &
826	1., &
827	tint, &
828	A8bd(1,1,loind), &
829	A8bd(1,1,hiind), &
830	A8cr)
831
832	! ... Put on model pressure levels
833
834	index = -9999
835	CALL pinterp (PLON, &
836	A8cr, &
837	ploc, &
838	nlevl, &
839	A8, &
840	pmid, &
841	PLEV, &
842	index, &
843	.false. )
844	! ...
845
846	CALL linintp (lonlev, &
847	0., &
848	1., &
849	tint, &
850	A9bd(1,1,loind), &
851	A9bd(1,1,hiind), &
852	A9cr)
853
854	! ... Put on model pressure levels
855
856	index = -9999
857	CALL pinterp (PLON, &
858	A9cr, &
859	ploc, &
860	nlevl, &
861	A9, &
862	pmid, &
863	PLEV, &
864	index, &
865	.false. )
866
867	END SUBROUTINE OZLIN_INTERP
868
869
870	#ifdef STRAT
871	SUBROUTINE SAD_READ(calday)
872	!----------------------------------------------------------------------
873	! ... Read sulfate data from CMIP5
874	! R. Valorso
875	!----------------------------------------------------------------------
876
877	USE SAD_COM
878	USE INCA_DIM
879	USE MOD_INCA_PARA
880	USE PRINT_INCA
881
882	IMPLICIT NONE
883
884	INCLUDE 'netcdf.inc'
885
886	!----------------------------------------------------------------------
887	! ... Dummy args
888	!----------------------------------------------------------------------
889	REAL, INTENT(IN) :: calday
890
891	!----------------------------------------------------------------------
892	! ... Local variables
893	!----------------------------------------------------------------------
894	INTEGER :: iret
895	INTEGER :: varid1,varid2,varid3,varid4,varid5
896	INTEGER :: oldlotime, oldhitime
897	INTEGER, DIMENSION(3) :: start, count
898	REAL :: PRES_bd_glo(nbp_glo,nlevs,2)
899	REAL :: SAD_bd_glo(nbp_glo,nlevs,2)
900	REAL :: MASS_bd_glo(nbp_glo,nlevs,2)
901	REAL :: VOL_bd_glo(nbp_glo,nlevs,2)
902	REAL :: RMEAN_bd_glo(nbp_glo,nlevs,2)
903
904	! ... Check to see if model time still bounded by data set
905	oldlotime = lotime
906	oldhitime = hitime
907	CALL FINDPLB (times, ntimes, calday, lotime)
908	IF ( cyclical ) THEN
909	hitime = MOD(lotime,ntimes) + 1
910	ELSE
911	hitime = lotime + 1
912	END IF
913
914	! ... Read new data if needed
915	IF ( hitime /= oldhitime ) THEN
916	loind = hiind
917	hiind = MOD( loind, 2) + 1
918
919	!$OMP MASTER
920	IF (is_mpi_root) THEN
921	iret = nf_inq_varid(ncids, 'PRES', varid1)
922	CALL check_err(iret, 'SAF_READ', 'problem to get id for pres')
923	iret = nf_inq_varid(ncids, 'SAD', varid2)
924	CALL check_err(iret, 'SAD_READ', 'problem to get id for sad')
925	iret = nf_inq_varid(ncids, 'MASS', varid3)
926	CALL check_err(iret, 'SAD_READ', 'problem to get id for mass')
927	iret = nf_inq_varid(ncids, 'VOL', varid4)
928	CALL check_err(iret, 'SAD_READ', 'problem to get id for vol')
929	iret = nf_inq_varid(ncids, 'RMEAN', varid5)
930	CALL check_err(iret, 'SAD_READ', 'problem to get id for rmean')
931
932	count(1) = nbp_glo
933	count(2) = nlevs
934	count(3) = 1
935	start(1) = 1
936	start(2) = 1
937
938	start(3) = hitime
939	iret = nf_get_vara_double(ncids, varid1, start, count, PRES_bd_glo(1,1,hiind))
940	CALL check_err(iret, 'SAD_READ', 'problem to read PRES_bd hiind')
941	iret = nf_get_vara_double(ncids, varid2, start, count, SAD_bd_glo(1,1,hiind))
942	CALL check_err(iret, 'SAD_READ', 'problem to read SAD_bd hiind')
943	iret = nf_get_vara_double(ncids, varid3, start, count, MASS_bd_glo(1,1,hiind))
944	CALL check_err(iret, 'SAD_READ', 'problem to read MASS_bd hiind')
945	iret = nf_get_vara_double(ncids, varid4, start, count, VOL_bd_glo(1,1,hiind))
946	CALL check_err(iret, 'SAD_READ', 'problem to read VOL_bd hiind')
947	iret = nf_get_vara_double(ncids, varid5, start, count, RMEAN_bd_glo(1,1,hiind))
948	CALL check_err(iret, 'SAD_READ', 'problem to read RMEAN_bd hiind')
949
950	ENDIF
951	!$OMP END MASTER
952	CALL scatter(PRES_bd_glo(:,:,hiind),PRES_bd(:,:,hiind))
953	CALL scatter(SAD_bd_glo(:,:,hiind),SAD_bd(:,:,hiind))
954	CALL scatter(MASS_bd_glo(:,:,hiind),MASS_bd(:,:,hiind))
955	CALL scatter(VOL_bd_glo(:,:,hiind),VOL_bd(:,:,hiind))
956	CALL scatter(RMEAN_bd_glo(:,:,hiind),RMEAN_bd(:,:,hiind))
957
958	IF ( lotime /= oldhitime ) THEN
959	!$OMP MASTER
960	IF (is_mpi_root) THEN
961	start(3) = lotime
962	iret = nf_get_vara_double(ncids, varid1, start, count, PRES_bd_glo(1,1,loind))
963	CALL check_err(iret, 'SAD_READ', 'problem to read PRES_bd loind')
964	iret = nf_get_vara_double(ncids, varid2, start, count, SAD_bd_glo(1,1,loind))
965	CALL check_err(iret, 'SAD_READ', 'problem to read SAD_bd loind')
966	iret = nf_get_vara_double(ncids, varid3, start, count, MASS_bd_glo(1,1,loind))
967	CALL check_err(iret, 'SAD_READ', 'problem to read MASS_bd loind')
968	iret = nf_get_vara_double(ncids, varid4, start, count, VOL_bd_glo(1,1,loind))
969	CALL check_err(iret, 'SAD_READ', 'problem to read VOL_bd loind')
970	iret = nf_get_vara_double(ncids, varid5, start, count, RMEAN_bd_glo(1,1,loind))
971	CALL check_err(iret, 'SAD_READ', 'problem to read RMEAN_bd loind')
972
973	ENDIF
974	!$OMP END MASTER
975	CALL scatter(PRES_bd_glo(:,:,loind),PRES_bd(:,:,loind))
976	CALL scatter(SAD_bd_glo(:,:,loind),SAD_bd(:,:,loind))
977	CALL scatter(MASS_bd_glo(:,:,loind),MASS_bd(:,:,loind))
978	CALL scatter(VOL_bd_glo(:,:,loind),VOL_bd(:,:,loind))
979	CALL scatter(RMEAN_bd_glo(:,:,loind),RMEAN_bd(:,:,loind))
980
981	END IF
982
983	END IF
984
985	END SUBROUTINE SAD_READ
986
987	SUBROUTINE SAD_INTERP(calday, pmid)
988	!----------------------------------------------------------------------
989	! ... Interpolate sulfate data from CMIP 5 on model time and on pressure levels
990	! R. Valorso
991	!----------------------------------------------------------------------
992
993	USE SAD_COM
994	USE INCA_DIM
995	IMPLICIT NONE
996
997	!----------------------------------------------------------------------
998	! ... Dummy args
999	!----------------------------------------------------------------------
1000	REAL, INTENT(IN) :: calday
1001	REAL, INTENT(IN) :: pmid(PLON,PLEV)
1002
1003	!----------------------------------------------------------------------
1004	! ... Local variables
1005	!----------------------------------------------------------------------
1006	REAL :: dt, dt1, tint
1007	REAL :: ploc(klons,nlevs)
1008	INTEGER :: lonlev
1009	INTEGER :: i
1010	INTEGER :: index(PLON,PLEV)
1011
1012	index = -9999
1013	lonlev = klons * nlevs
1014
1015	!----------------------------------------------------------------------
1016	! ... First interpolate linearly on current model time
1017
1018	IF ( times(hitime) < times(lotime) ) THEN
1019	dt = 365. - times(lotime) + times(hitime)
1020	IF (calday <= times(hitime) ) THEN
1021	dt1 = 365. - times(lotime) + calday
1022	ELSE
1023	dt1 = calday - times(lotime)
1024	END IF
1025	ELSE
1026	dt = times(hitime) - times(lotime)
1027	dt1 = calday - times(lotime)
1028	END IF
1029
1030	tint = dt1/dt
1031
1032	CALL linintp (lonlev, &
1033	0., &
1034	1., &
1035	tint, &
1036	PRES_bd(1,1,loind), &
1037	PRES_bd(1,1,hiind), &
1038	PRES_cr)
1039
1040	DO i = 1, klons
1041	ploc(i,:) = PRES_cr(i,:)
1042	ENDDO
1043
1044	! ...
1045
1046	CALL linintp (lonlev, &
1047	0., &
1048	1., &
1049	tint, &
1050	SAD_bd(1,1,loind), &
1051	SAD_bd(1,1,hiind), &
1052	SAD_cr)
1053
1054	! ... Put on model pressure levels
1055
1056	index = -9999
1057	CALL pinterp (PLON, &
1058	SAD_cr, &
1059	ploc, &
1060	nlevs, &
1061	SAD, &
1062	pmid, &
1063	PLEV, &
1064	index, &
1065	.false. )
1066	! ...
1067
1068	CALL linintp (lonlev, &
1069	0., &
1070	1., &
1071	tint, &
1072	MASS_bd(1,1,loind), &
1073	MASS_bd(1,1,hiind), &
1074	MASS_cr)
1075
1076	! ... Put on model pressure levels
1077
1078	index = -9999
1079	CALL pinterp (PLON, &
1080	MASS_cr, &
1081	ploc, &
1082	nlevs, &
1083	MASS, &
1084	pmid, &
1085	PLEV, &
1086	index, &
1087	.false. )
1088
1089	! ...
1090
1091	CALL linintp (lonlev, &
1092	0., &
1093	1., &
1094	tint, &
1095	VOL_bd(1,1,loind), &
1096	VOL_bd(1,1,hiind), &
1097	VOL_cr)
1098
1099	! ... Put on model pressure levels
1100
1101	index = -9999
1102	CALL pinterp (PLON, &
1103	VOL_cr, &
1104	ploc, &
1105	nlevs, &
1106	VOL, &
1107	pmid, &
1108	PLEV, &
1109	index, &
1110	.false. )
1111	! ...
1112
1113	CALL linintp (lonlev, &
1114	0., &
1115	1., &
1116	tint, &
1117	RMEAN_bd(1,1,loind), &
1118	RMEAN_bd(1,1,hiind), &
1119	RMEAN_cr)
1120
1121	! ... Put on model pressure levels
1122
1123	index = -9999
1124	CALL pinterp (PLON, &
1125	RMEAN_cr, &
1126	ploc, &
1127	nlevs, &
1128	RMEAN, &
1129	pmid, &
1130	PLEV, &
1131	index, &
1132	.false. )
1133
1134
1135	END SUBROUTINE SAD_INTERP
1136	#endif
1137	#endif
1138	#if !defined(GES) && !defined(DUSS)
1139	SUBROUTINE BOUNDSPC (dtinv, mmr, pmid, temp)
1140	!----------------------------------------------------------------------
1141	! ... Impose boundary conditions for chemical tracers
1142	! Didier Hauglustaine, IPSL, 1999.
1143	!----------------------------------------------------------------------
1144
1145	USE O3CLIM_COM
1146	USE SPECIES_NAMES
1147	USE OXYDANT_COM
1148	USE AIRPLANE_SRC, ONLY : itrop
1149	USE CHEM_MODS, ONLY : adv_mass, invariants
1150	USE INCA_DIM
1151	#ifdef STRAT
1152	USE LGLIVED_MOD
1153	#endif
1154	USE PARAM_CHEM
1155	USE RATE_INDEX_MOD
1156
1157	IMPLICIT NONE
1158
1159	!----------------------------------------------------------------------
1160	! ... Dummy args
1161	!----------------------------------------------------------------------
1162	REAL, INTENT(IN) :: dtinv
1163	REAL, INTENT(IN) :: pmid(PLON,PLEV)
1164	REAL, INTENT(IN) :: temp(PLON,PLEV)
1165	REAL, INTENT(INOUT) :: mmr(PLON,PLEV,PCNST)
1166
1167	!----------------------------------------------------------------------
1168	! ... Local variables
1169	!----------------------------------------------------------------------
1170	INTEGER :: i, k
1171	INTEGER :: ktrop, krelax, k380
1172	REAL :: taurelax, facrelax
1173	REAL :: delt
1174	REAL, PARAMETER :: dry_mass = 28.966
1175
1176	REAL :: theta(PLON,PLEV)
1177
1178	delt = 1./dtinv
1179	taurelax = 86400.*10. ! (10 days)
1180	facrelax = 1. - EXP( -delt / taurelax )
1181	theta(:,:) = temp(:,:)(1.e5/pmid(:,:))*0.286
1182
1183	DO i = 1, PLON
1184
1185	ktrop = itrop(i)
1186
1187	DO k = ktrop, PLEV
1188	IF (theta(i,k) >= 380.) THEN
1189	k380 = k
1190	EXIT
1191	ENDIF
1192	ENDDO
1193
1194	! krelax = MIN(ktrop+2,PLEV) ! 2 levels above tropopause
1195	krelax = k380 ! 380K theta surface
1196
1197
1198	! ... Impose NOy to climatologies
1199	! mmr(i,krelax:PLEV,id_NO) = nomzrt(i,krelax:PLEV) * adv_mass(id_NO)/dry_mass
1200	#ifdef AERONLY
1201	mmr(i,krelax:PLEV,id_NO2) = invariants(i,krelax:PLEV, inv_NO2INV)/invariants(i,krelax:PLEV, inv_M) * adv_mass(id_NO2)/dry_mass
1202	mmr(i,krelax:PLEV,id_HNO3) = invariants(i,krelax:PLEV, inv_HNO3INV)/invariants(i,krelax:PLEV, inv_M) * adv_mass(id_HNO3)/dry_mass
1203	#else
1204	! ... Relax O3 to climatologies
1205	if (trim(flag_o3) .eq. "o3clim") then
1206	mmr(i,krelax:PLEV,id_O3) = o3clim(i,krelax:PLEV)
1207	elseif (trim(flag_o3) .eq. "o3lin") then
1208	mmr(i,krelax:PLEV,id_O3) = mmr(i,ktrop:PLEV,id_O3L)
1209	endif
1210
1211	! ... Impose inert O3 to O3 above tropopause
1212	mmr(i,ktrop:PLEV,id_O3I) = mmr(i,ktrop:PLEV,id_O3)
1213
1214	! ... Impose stratospheric O3 to O3 above tropopause
1215	mmr(i,ktrop:PLEV,id_O3S) = mmr(i,ktrop:PLEV,id_O3)
1216
1217	#endif
1218	! ... Impose CH4
1219	! mmr(i,:,id_CH4) = 1760E-9 * 16./dry_mass ! forcage CH4 pour IPCC2000
1220	! mmr(i,:,id_CH4) = 791.6E-9 * 16./dry_mass ! forcage CH4 pour IPCC2000
1221
1222	#ifdef STRAT
1223	! Stratospheric Version
1224	! Impose long_lived species at the surface CMIP5
1225
1226	mmr(i,1,id_CH4) = conc_cfc(1) * adv_mass(id_CH4) / dry_mass
1227	mmr(i,1,id_N2O) = conc_cfc(2) * adv_mass(id_N2O) / dry_mass
1228
1229	mmr(i,1,id_CFC11) = conc_cfc(3) * adv_mass(id_CFC11) / dry_mass
1230	mmr(i,1,id_CFC12) = conc_cfc(4) * adv_mass(id_CFC12) / dry_mass
1231	mmr(i,1,id_CFC113) = conc_cfc(5) * adv_mass(id_CFC113) / dry_mass
1232	mmr(i,1,id_HCFC22) = conc_cfc(6) * adv_mass(id_HCFC22) / dry_mass
1233	mmr(i,1,id_HFC134a) = conc_cfc(7) * adv_mass(id_HFC134a) / dry_mass
1234	mmr(i,1,id_HCFC141) = conc_cfc(8) * adv_mass(id_HCFC141) / dry_mass
1235	mmr(i,1,id_HCFC142) = conc_cfc(9) * adv_mass(id_HCFC142) / dry_mass
1236	mmr(i,1,id_Ha1301) = conc_cfc(10) * adv_mass(id_Ha1301) / dry_mass
1237	mmr(i,1,id_Ha1211) = conc_cfc(11) * adv_mass(id_Ha1211) / dry_mass
1238	mmr(i,1,id_MCF) = conc_cfc(12) * adv_mass(id_MCF) / dry_mass
1239	mmr(i,1,id_CCl4) = conc_cfc(13) * adv_mass(id_CCl4) / dry_mass
1240	mmr(i,1,id_CH3Cl) = conc_cfc(14) * adv_mass(id_CH3Cl) / dry_mass
1241	mmr(i,1,id_CH3Br) = conc_cfc(15) * adv_mass(id_CH3Br) / dry_mass
1242	#endif
1243
1244	END DO
1245
1246	END SUBROUTINE BOUNDSPC
1247	#endif
1248	SUBROUTINE FDTROPOPAUSE (nx, ny, nz, pmid, temp)
1249	!****************************************************************************
1250	! Purpose
1251	! -------
1252	!
1253	! Calculates 2D fields of thermal tropopause pressures and tropopause
1254	! level indices for given 3D fields of temperature and pressure.
1255	!
1256	! Methods
1257	! -------
1258	! - Subroutine stattrop calculates the thermal tropopause pressure (Pa)
1259	! for a 1D column (sounding) of pressures (Pa) and temperatures (K).
1260	!
1261	! - The program tropo_test also corrects for non-sense occuring at high
1262	! latitudes > 60N (60S).
1263	! Theese problems occur at high latitudes on the winter hemisphere because
1264	! of the cold stratosphere and the consequently weak troposphere-
1265	! stratosphere transition of the thermal lapse rate -dT/dz which is
1266	! used to determine the thermal tropopause.
1267	! The problem is much less severe on the winter hemisphere.
1268	!
1269	! Here a fix is implemented applying a factor which scales the
1270	! tropopause pressures at latitudes LAT > 60N (60S) to the
1271	! mean tropopause pressure at 60N (60S), whenever the mean
1272	! tropopause pressure at LAT is lower than that at 60N (60S).
1273	! This idea was formulated in the ECHAM4/CHEM routine xttphwmo by
1274	! D. Nodorp, Ch. Land, B. Steil and R. Hein.
1275	!
1276	! Programmed by Dominik Brunner, ETHZ, Switzerland V.01, 10 Aug 2000
1277	! Modified by Didier Hauglustaine, IPSL, for LMDZ/INCA, Oct 2000.
1278	!
1279	!*****************************************************************************
1280
1281	USE AIRPLANE_SRC, ONLY : itrop, ptrop
1282	USE INCA_DIM
1283	USE MOD_INCA_PARA, ONLY : ij_begin,ij_end, &
1284	nbp_para_begin,nbp_para_end,mpi_rank, &
1285	plon_omp_begin, plon_omp_end
1286	USE MOD_GRID_INCA
1287	USE MOD_INCA_OMP_DATA
1288
1289	IMPLICIT NONE
1290	!----------------------------------------------------------------------
1291	! ... Dummy args
1292	!----------------------------------------------------------------------
1293	REAL, INTENT(IN) :: pmid(PLON_GLO,PLEV)
1294	REAL, INTENT(IN) :: temp(PLON_GLO,PLEV)
1295	INTEGER, INTENT(IN) :: nx, ny, nz
1296
1297	!----------------------------------------------------------------------
1298	! ... Local variables
1299	!----------------------------------------------------------------------
1300	REAL,DIMENSION(nx,ny,nz) :: p3 ! pressures at model layers (full levels) (Pa)
1301	REAL,DIMENSION(nx,ny,nz) :: t3 ! temp. at model layer (full levels) (K)
1302	INTEGER :: i,j,k,l ! lon, lat and vertical indices
1303	REAL, DIMENSION(nx,ny) :: ptropd ! 2D field of tropopause pressures
1304	INTEGER, DIMENSION(nx,ny) :: itropd ! 2D field of tropopause level indices
1305	REAL, DIMENSION(nx,ny) :: itropdr ! 2D field of tropopause level indices
1306	REAL, DIMENSION(PLON_GLO) :: itropr
1307	REAL, DIMENSION(nz) :: p1,t1 ! 1D columns of pressures & temperatures
1308	INTEGER :: jlat60S, jlat60N ! latitude indices for 60N and 60S
1309	REAL, PARAMETER :: rinval=-999.00 ! invalid real number
1310	INTEGER :: l300hPa ! index of layer containing 300 hPa level
1311	INTEGER :: l120hPa ! index of layer containing 120 hPa level
1312	REAL :: ptpm60,ptpmlat,ptpsum ! mean and integrated TP pressures
1313	REAL :: rscale ! scaling factor to correct non-sense
1314	REAL :: deltalat
1315	REAL :: ptrop_glo(PLON_GLO)
1316	INTEGER :: itrop_glo(PLON_GLO)
1317	INTEGER :: nbbeg_loc,nbend_loc
1318
1319
1320
1321	deltalat = 180./FLOAT(ny-1)
1322
1323	!T and p to 2D grid
1324	CALL grid1dto2d_glo(pmid,p3)
1325	CALL grid1dTo2d_glo(temp,t3)
1326
1327	!indices of lat=60N and 60S
1328	jlat60N=FLOOR(30./deltalat)+1
1329	jlat60S=ny-jlat60N+1
1330
1331	!Calc index of model layer (typically) containing the 300 hPa level
1332	!loop from surface to top
1333	DO l=1,nz-1
1334	IF (0.5*(LOG(p3(nx/2,ny/2,l))+LOG(p3(nx/2,ny/2,l+1))).LT.LOG(30000.)) &
1335	GOTO 100
1336	ENDDO
1337	100 CONTINUE
1338	l300hPa=l
1339
1340	!Calc index of model layer (typically) containing the 120 hPa level
1341	!loop from surface to top
1342	DO l=1,nz-1
1343	IF (0.5*(LOG(p3(nx/2,ny/2,l))+LOG(p3(nx/2,ny/2,l+1))).LT.LOG(12000.)) &
1344	GOTO 101
1345	ENDDO
1346	101 CONTINUE
1347	l120hPa=l
1348
1349	!loop over southern hemisphere and get tropopause (loop starts at equator!!)
1350	DO j=ny/2,ny
1351	ptpsum=0. ! initialize sum of tropopause pressures
1352	DO i=1,nx
1353	!reverse order of levels because stattrop
1354	!requires first index at top of atmosphere
1355	DO l=1,nz
1356	p1(l)=p3(i,j,nz+1-l)
1357	t1(l)=t3(i,j,nz+1-l)
1358	ENDDO
1359	CALL stattrop(p1,t1,nz,ptropd(i,j),itropd(i,j))
1360	itropd(i,j)=nz-itropd(i,j)+1 ! reverse tropopause level index back
1361
1362	!if no tropopause was found (happens only rarely)
1363	IF (ptropd(i,j).EQ.rinval) THEN
1364	IF (j.GT.jlat60S) THEN
1365	ptropd(i,j)=30000. ! set to 300 hPa at high latitudes
1366	itropd(i,j)=l300hPa ! index of layer containing 300 hPa level
1367	ELSE
1368	ptropd(i,j)=12000. ! set to 120 hPa at low latitudes
1369	itropd(i,j)=l120hPa ! index of layer containing 120 hPa level
1370	END IF
1371	ENDIF
1372
1373	!sum up tropopause levels to calculate mean
1374	ptpsum=ptpsum+ptropd(i,j)
1375	ENDDO
1376
1377	!Correct nonsense at high latitudes (correction is only
1378	!necessary during southern hemisphere winter from Apr until Sep)
1379	IF (j.EQ.jlat60S) ptpm60=ptpsum/float(nx)
1380
1381	IF (j.GT.jlat60S) THEN
1382	ptpmlat=ptpsum/float(nx)
1383	!correct ptropd if mean tropopause pressure at this latitude
1384	!is lower than at 60S
1385	IF (ptpmlat.LT.ptpm60) THEN
1386	rscale=ptpm60/ptpmlat
1387	DO i=1,nx
1388	ptropd(i,j)=ptropd(i,j)*rscale
1389	!find corresponding level indices (start loop at surface)
1390	DO l=1,nz-1
1391	IF (0.5*(LOG(p3(i,j,l))+LOG(p3(i,j,l+1))) &
1392	.LT.LOG(ptropd(i,j))) GOTO 102
1393	ENDDO
1394	102 itropd(i,j)=l
1395	ENDDO
1396	ENDIF
1397	ENDIF
1398	ENDDO
1399
1400	!loop over northern hemisphere and get tropopause (loop starts at equator!!)
1401	DO j=ny/2,1,-1
1402	ptpsum=0. ! initialize sum of tropopause pressures
1403	DO i=1,nx
1404	!reverse order of levels because stattrop
1405	!requires first index at top of atmosphere
1406	DO l=1,nz
1407	p1(l)=p3(i,j,nz+1-l)
1408	t1(l)=t3(i,j,nz+1-l)
1409	ENDDO
1410	CALL stattrop(p1,t1,nz,ptropd(i,j),itropd(i,j))
1411	itropd(i,j)=nz-itropd(i,j)+1 ! reverse tropopause level index back
1412
1413	!if no tropopause was found (happens only rarely)
1414	IF (ptropd(i,j).EQ.rinval) THEN
1415	IF (j.LT.jlat60N) THEN
1416	ptropd(i,j)=30000. ! set to 300 hPa for high latitudes
1417	itropd(i,j)=l300hPa ! index of layer containing 300 hPa level
1418	ELSE
1419	ptropd(i,j)=12000. ! set to 120 hPa for low latitudes
1420	itropd(i,j)=l120hPa ! index of layer containing 120 hPa level
1421	END IF
1422	ENDIF
1423
1424	!sum up tropopause levels to calculate mean
1425	ptpsum=ptpsum+ptropd(i,j)
1426	ENDDO
1427
1428	!Correct nonsense at high latitudes (correction is only
1429	!necessary during nouthern hemisphere winter from Oct-Mar)
1430	IF (j.EQ.jlat60N) ptpm60=ptpsum/float(nx)
1431
1432	IF (j.LT.jlat60N) THEN
1433	ptpmlat=ptpsum/float(nx)
1434	!correct ptropd if mean tropopause pressure at this latitude
1435	!is lower than at 60S
1436	IF (ptpmlat.LT.ptpm60) THEN
1437	rscale=ptpm60/ptpmlat
1438	DO i=1,nx
1439	ptropd(i,j)=ptropd(i,j)*rscale
1440	!find corresponding level indices (start loop at surface)
1441	DO l=1,nz-1
1442	IF (0.5*(LOG(p3(i,j,l))+LOG(p3(i,j,l+1))) &
1443	.LT.LOG(ptropd(i,j))) GOTO 103
1444	ENDDO
1445	103 itropd(i,j)=l
1446	ENDDO
1447	ENDIF
1448	ENDIF
1449	ENDDO
1450
1451	! borne pour le passage GLO -> LOC
1452	nbbeg_loc = nbp_para_begin(mpi_rank)+plon_omp_begin-1
1453	nbend_loc = nbp_para_begin(mpi_rank)+plon_omp_end-1
1454	!Back to physiq grid 1D
1455	CALL grid2dto1d_glo(ptropd,ptrop_glo)
1456	ptrop=ptrop_glo(nbbeg_loc:nbend_loc)
1457
1458
1459	itropdr = FLOAT(itropd)
1460	CALL grid2dto1d_glo(itropdr,itropr)
1461	itrop_glo = INT(itropr)
1462
1463	itrop=itrop_glo(nbbeg_loc:nbend_loc)
1464
1465	END SUBROUTINE FDTROPOPAUSE
1466
1467	SUBROUTINE PINTERP(ncol, &
1468	fieldin, &
1469	pin, &
1470	nin, &
1471	fieldout, &
1472	pout, &
1473	nout, &
1474	index, &
1475	entered)
1476	!-----------------------------------------------------------------------
1477	! ... Interpolates on model pressure levels
1478	! Stacy Walters and Didier Hauglustaine, 1999.
1479	!-----------------------------------------------------------------------
1480
1481	IMPLICIT NONE
1482
1483	!-----------------------------------------------------------------------
1484	! ... Dummy args
1485	!-----------------------------------------------------------------------
1486	LOGICAL, INTENT(in) :: entered
1487	INTEGER, INTENT(in) :: nin
1488	INTEGER, INTENT(in) :: nout
1489	INTEGER, INTENT(in) :: ncol
1490	INTEGER, INTENT(inout) :: INDEX(ncol,nout)
1491	REAL, INTENT(in) :: pin(ncol,nin)
1492	REAL, INTENT(in) :: fieldin(ncol,nin)
1493	REAL, INTENT(in) :: pout(ncol,nout)
1494	REAL, INTENT(out) :: fieldout(ncol,nout)
1495
1496	!-----------------------------------------------------------------------
1497	! ... Local variables
1498	!-----------------------------------------------------------------------
1499	INTEGER :: i, j, k
1500	REAL :: deltp(ncol,nout)
1501	REAL :: pinratio(ncol,nin)
1502
1503	IF ( .NOT. entered) THEN
1504
1505	DO i = 1, ncol
1506	DO k = 1, nout
1507
1508	DO j = nin-1, 1, -1
1509	IF (pout(i,k) <= pin(i,j)) THEN
1510	IF (pout(i,k) > pin(i,j+1)) THEN
1511	index(i,k) = j
1512	ENDIF
1513	ENDIF
1514	ENDDO
1515
1516	ENDDO
1517	ENDDO
1518
1519	ENDIF
1520
1521	DO i = 1, ncol
1522	DO j = 1, nin-1
1523	pinratio(i,j) = 1./LOG(pin(i,j)/pin(i,j+1))
1524	END DO
1525	END DO
1526
1527
1528	DO i = 1, ncol
1529	DO k = 1, nout
1530
1531
1532	IF (pout(i,k) <= pin(i,nin)) THEN
1533	fieldout(i,k) = fieldin(i,nin)
1534	ELSE IF (pout(i,k) >= pin(i,1)) THEN
1535	fieldout(i,k) = fieldin(i,1)
1536	ELSE
1537
1538	deltp(i,k) = LOG(pout(i,k)/pin(i,index(i,k)+1)) * pinratio(i,index(i,k))
1539
1540	fieldout(i,k) = fieldin(i,index(i,k)+1) &
1541	+ deltp(i,k) * (fieldin(i,index(i,k))-fieldin(i,index(i,k)+1))
1542
1543	ENDIF
1544
1545	ENDDO
1546	ENDDO
1547
1548	END SUBROUTINE PINTERP
1549
1550	SUBROUTINE stattrop(pfull, tfull, nlev, ptropd, itropd)
1551
1552	!**********************************************************************************
1553	!
1554	! input : p, t real(nlev)
1555	! input : nlev real
1556	! output: ptropd real
1557	! output: itropd integer
1558	!
1559	! programmed by Dominik Brunner V1.0 Aug 2000
1560	!
1561	! built upon routine stattrop by Peter van Velthoven,
1562	! KNMI, The Netherlands
1563	! and on the ECHAM4/CHEM routine xttphwmo by
1564	! Thomas Reichle, Christine Land, B. Steil and R. Hein, DLR
1565	!
1566	! purpose
1567	! -------
1568	! stattrop computes the pressure (Pa) at the thermal (static) tropopause
1569	! (TP) for a 1D column (sounding) of pressures and temperatures following
1570	! the definition of the height of the tropopause as postulated by WMO (1957).
1571	!
1572	! ATTENTION: In the current formulation of the program the first
1573	! level (index 1) must be at the top of the atmosphere and the
1574	! last level (index nlev) at the surface
1575	!
1576	! interface
1577	! ---------
1578	! call stattrop(pfull, tfull, nlev, ptropd, itropd)
1579	! - Input
1580	! nlev : number of vertical levels
1581	! pfull: pressure in 1D column at nlev full levels (layers)
1582	! tfull: temperature in 1D column at nlev full levels
1583	! - Output
1584	! ptropd: height of the tropopause in Pa
1585	! itropd: index of layer containing the tropopause
1586	!
1587	! methods
1588	! -------
1589	! - Lapse rate gamma = -dT/dz
1590	! Using the hydrostatic approximation
1591	!
1592	! dz = -dp/(grho) = -dp/p R/g * T = -dlnp * R/g * T
1593	!
1594	! we get -dT/dz = dT/T * g/R *1/dlnp = dlnT/dlnp
1595	!
1596	! - Variables are assumed to vary linearly in log(pressure)
1597	! - The tropopause is the lowest level above 450 hPa fullfilling
1598	! the WMO criterium. If ptropd is < 85 hPa it is set to 85 hPa.
1599	! If no tropopause is found ptropd is set to -999.
1600	!
1601	! references
1602	! ----------
1603	!
1604	! - WMO (1992): International meteorological vocabulary, Genf, 784pp.:
1605	!
1606	! 1. The first tropopause is defined as the lowest level at which
1607	! the lapse rate decreases to 2 deg C per kilometer or less,
1608	! provided also the average lapse rate between this level and
1609	! all higher levels within 2 kilometers does not exceed 2 deg C
1610	!
1611	! - Randel WJ, Wu F, Gaffen DJ, Interannual variability of the tropical
1612	! tropopause derived from radiosonde data and NCEP reanalyses,
1613	! JOURNAL OF GEOPHYSICAL RESEARCH, 105, 15509-15523, 2000.
1614	!
1615	! The following webpage clearifies the calculation of the tropopause
1616	! in the NCEP reanalysis: http://dss.ucar.edu/pub/reanalysis/FAQ.html
1617	!
1618	! For comparison NCEP reanalysis tropopause pressures can be obtained
1619	! from http://www.cdc.noaa.gov/cdc/reanalysis/reanalysis.shtml
1620	!
1621	!**********************************************************************************
1622	USE INCA_DIM
1623
1624	IMPLICIT NONE
1625	INTEGER nlev
1626	REAL pfull(nlev), tfull(nlev)
1627	REAL ptropd, p2km, lapseavg, lapsesum
1628	INTEGER ilev, ilev1, itropd, count, l
1629
1630	REAL lapse(PLEV), phalf(PLEV)
1631	REAL grav, rgas, const, lapsec, rinval
1632	PARAMETER (grav=9.80665, rgas=287.04, &
1633	const=1000.*grav/rgas, &
1634	lapsec=2.0, rinval=-999.0 )
1635	! min. and max. allowed tropopause pressure
1636	REAL pmin, pmax
1637	PARAMETER (pmin=8500.,pmax=45000.)
1638	! deltaz is layer thickness used for second tropopause criterium
1639	REAL deltaz
1640	PARAMETER (deltaz=2.)
1641	LOGICAL found
1642	REAL p1, p2 , weight
1643
1644	ptropd = rinval
1645	itropd = 0
1646	found = .FALSE.
1647
1648	! Loop from model top to surface and calculate lapse rate gamma=-dT/dz (K/km)
1649	DO ilev = 1, nlev-1
1650	ilev1 = ilev + 1
1651	lapse(ilev) = LOG(tfull(ilev))- LOG(tfull(ilev1))
1652	lapse(ilev) = lapse(ilev) / &
1653	(LOG(pfull(ilev))- LOG(pfull(ilev1)))
1654	lapse(ilev) = const * lapse(ilev)
1655	phalf(ilev) = (pfull(ilev) + pfull(ilev1))*0.5
1656	END DO
1657	! Loop from surface to top to find (lowest) tropopause
1658	DO ilev = nlev-2, 1, -1
1659	! ** 1st test: -dT/dz less than 2 K/km and pressure LT pmax? **
1660	!Modified 07/2001 -- D. Hauglustaine
1661	IF (lapse(ilev).LT.lapsec.AND.pfull(ilev).LT.pmax.AND.ilev.GE.5) THEN
1662	IF (.NOT.found) THEN
1663	! Interpolate tropopause pressure linearly in log(p)
1664	p1 = LOG(phalf(ilev))
1665	p2 = LOG(phalf(ilev+1))
1666	!Modified 09/2001 -- L. Jourdain
1667	IF ( (lapse(ilev).NE.lapse(ilev+1)) .AND. (lapse(ilev+1).GE.lapsec) ) THEN
1668	weight = (lapsec - lapse(ilev)) / &
1669	(lapse(ilev+1)-lapse(ilev))
1670	ptropd = EXP(p1 + weight * (p2-p1)) ! tropopause pressure
1671	ELSE
1672	ptropd = phalf(ilev)
1673	END IF
1674	! *** 2nd test: average -dT/dz in layer deltaz above TP must not be greater than 2K/km
1675	p2km = EXP(LOG(ptropd)-deltaz*const/tfull(ilev)) ! press 2 km above TP
1676	lapsesum = 0.0 ! init avg. lapse rate 2 km above TP
1677	lapseavg=1.e9
1678	count = 0
1679	DO l=ilev,1,-1
1680	IF (phalf(l).LT.p2km) GOTO 100
1681	lapsesum=lapsesum+lapse(l)
1682	count=count+1
1683	lapseavg=lapsesum/float(count)
1684	END DO
1685	100 CONTINUE
1686	found=lapseavg.LT.lapsec
1687	! Discard tropopause?
1688	IF (.NOT.found) THEN
1689	ptropd=rinval
1690	ELSE
1691	ptropd=MAX(ptropd,pmin) ! ptropd must be GE 85 hPa
1692	itropd=ilev ! assign level index
1693	END IF
1694	END IF
1695	END IF
1696	END DO
1697	RETURN
1698	END SUBROUTINE stattrop
1699
1700
1701
1702	subroutine CALC_PV(latitude_deg,paprs,pplay,ts,t,rot)
1703	! This subroutine consists in deriving Ertel's potential vorticity (PV) in
1704	! potential vorticity units (PVU), i.e. 1 PVU = 1E-6 K.m^2/kg/s.
1705	! It accounts for the vertical potential temperature gradient (=static stability)
1706	! and the (almost) horizontal wind vorticity.
1707	! The equation and clear explanations are available in Gettelman et al. (2011, Rev. Geophys),
1708	! a review paper on the extratropical UTLS.
1709	! Added by Yann Cohen, November 2019. ! Created by Yann Cohen
1710
1711	USE CONST_MOD
1712	USE XIOS_INCA
1713	IMPLICIT NONE
1714
1715	!-------------------------------------------------------------------------------
1716	! Arguments:
1717	REAL, INTENT(IN) :: latitude_deg(PLON) ! latitude
1718	REAL, INTENT(IN) :: paprs(PLON,PLEVP) !--- Cells-edges pressure
1719	REAL, INTENT(IN) :: pplay(PLON,PLEV) !--- Cells-centers pressure
1720	REAL, INTENT(IN) :: ts(PLON) !--- Surface temperature
1721	REAL, INTENT(IN) :: t(PLON,PLEV) !--- Cells-centers temperature
1722	REAL, INTENT(IN) :: rot(PLON,PLEV) !--- Cells-centers relative vorticity
1723	!-------------------------------------------------------------------------------
1724	! Local variables:
1725	include "YOMCST_I.h"
1726	INTEGER :: i, k
1727	REAL :: al,al2top
1728	REAL,PARAMETER :: preff=101325 ! Unit: Pa.
1729	REAL, DIMENSION(PLON,PLEVP):: tpot_edg
1730	REAL, DIMENSION(PLON,PLEV) :: tpot_cen
1731	REAL, DIMENSION(PLON,PLEV) :: PV
1732	!-------------------------------------------------------------------------------
1733	PV(:,:) = 0.
1734
1735	!--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES
1736	DO i = 1,PLON
1737	tpot_cen(i,1) = t(i,1)(preff/pplay(i,1))*RKAPPA
1738	tpot_edg(i,1) = ts(i) (preff/paprs(i,1))*RKAPPA
1739
1740	DO k=2,PLEV
1741
1742	tpot_cen(i,k) = t(i,k)(preff/pplay(i,k))*RKAPPA
1743	al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
1744	tpot_edg(i,k) = (t(i,k-1)+al(t(i,k)-t(i,k-1)))(preff/paprs(i,k))**RKAPPA
1745
1746	END DO
1747
1748	! Last ingredient for the PV field on the whole vertical grid:
1749	! linear extrapolation up to the summit, for cell-edges potential temperature.
1750	! The subsequent Theta vertical gradient is then used for the PV calculation at the top (full) level.
1751	! al2top=LOG(paprs(i,PLEV)/paprs(i,PLEVP))/LOG(paprs(i,PLEV)/pplay(i,PLEV))
1752	! tpot_edg(i,PLEVP) = tpot_edg(i,PLEV) + al2top*(tpot_cen(i,PLEV)-tpot_edg(i,PLEV))
1753
1754	END DO
1755
1756
1757	!--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
1758	DO i = 1, PLON
1759	DO k= 1, PLEV-1
1760	PV(i,k)=-1.E6RG(rot(i,k)+2.ROMEGASIN(latitude_deg(i)*RPI/180.))&
1761	* (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k))
1762	END DO
1763	END DO
1764
1765
1766	CALL xios_inca_change_context("inca")
1767	CALL xios_inca_send_field("PV", PV)
1768	CALL xios_inca_change_context("LMDZ")
1769
1770	end subroutine CALC_PV

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

Context Navigation

source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/src/INCA_SRC/set_ub_vals.F90 @ 6610

Download in other formats: