Context Navigation

calfis.f90 @ 5182

Last change on this file since 5182 was 5182, checked in by abarral, 10 days ago
(WIP) Replace REPROBUS CPP KEY by logical properly name modules
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:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 20.3 KB

Line
1	! $Id: calfis.f90 5182 2024-09-10 14:25:29Z abarral $
2
3
4
5	SUBROUTINE calfis(lafin, &
6	jD_cur, jH_cur, &
7	pucov, &
8	pvcov, &
9	pteta, &
10	pq, &
11	pmasse, &
12	pps, &
13	pp, &
14	ppk, &
15	pphis, &
16	pphi, &
17	pducov, &
18	pdvcov, &
19	pdteta, &
20	pdq, &
21	flxw, &
22	pdufi, &
23	pdvfi, &
24	pdhfi, &
25	pdqfi, &
26	pdpsfi)
27
28	! Auteur : P. Le Van, F. Hourdin
29	! .........
30	USE lmdz_infotrac, ONLY: nqtot, tracers
31	USE control_mod, ONLY: planet_type, nsplit_phys
32	USE callphysiq_mod, ONLY: call_physiq
33	USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
34	USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
35	USE comvert_mod, ONLY: preff, presnivs
36	USE lmdz_iniprint, ONLY: lunout, prt_level
37	USE lmdz_ssum_scopy, ONLY: scopy, ssum
38	USE lmdz_comgeom2
39	USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
40	USE lmdz_paramet
41
42	IMPLICIT NONE
43	!=======================================================================
44
45	! 1. rearrangement des tableaux et transformation
46	! variables dynamiques > variables physiques
47	! 2. calcul des termes physiques
48	! 3. retransformation des tendances physiques en tendances dynamiques
49
50	! remarques:
51	! ----------
52
53	! - les vents sont donnes dans la physique par leurs composantes
54	! naturelles.
55	! - la variable thermodynamique de la physique est une variable
56	! intensive : T
57	! pour la dynamique on prend T * ( preff / p(l) ) **kappa
58	! - les deux seules variables dependant de la geometrie necessaires
59	! pour la physique sont la latitude pour le rayonnement et
60	! l'aire de la maille quand on veut integrer une grandeur
61	! horizontalement.
62	! - les points de la physique sont les points scalaires de la
63	! la dynamique; numerotation:
64	! 1 pour le pole nord
65	! (jjm-1)*iim pour l'interieur du domaine
66	! ngridmx pour le pole sud
67	! ---> ngridmx=2+(jjm-1)*iim
68
69	! Input :
70	! -------
71	! pucov covariant zonal velocity
72	! pvcov covariant meridional velocity
73	! pteta potential temperature
74	! pps surface pressure
75	! pmasse masse d'air dans chaque maille
76	! pts surface temperature (K)
77	! callrad clef d'appel au rayonnement
78
79	! Output :
80	! --------
81	! pdufi tendency for the natural zonal velocity (ms-1)
82	! pdvfi tendency for the natural meridional velocity
83	! pdhfi tendency for the potential temperature
84	! pdtsfi tendency for the surface temperature
85
86	! pdtrad radiative tendencies \ both input
87	! pfluxrad radiative fluxes / and output
88
89	!=======================================================================
90
91	!-----------------------------------------------------------------------
92
93	! 0. Declarations :
94	! ------------------
95
96	INTEGER :: ngridmx
97	PARAMETER(ngridmx = 2 + (jjm - 1) * iim - 1 / jjm)
98
99	! Arguments :
100	! -----------
101	LOGICAL, INTENT(IN) :: lafin ! .TRUE. for the very last CALL to physics
102	REAL, INTENT(IN) :: jD_cur, jH_cur
103	REAL, INTENT(IN) :: pvcov(iip1, jjm, llm) ! covariant meridional velocity
104	REAL, INTENT(IN) :: pucov(iip1, jjp1, llm) ! covariant zonal velocity
105	REAL, INTENT(IN) :: pteta(iip1, jjp1, llm) ! potential temperature
106	REAL, INTENT(IN) :: pmasse(iip1, jjp1, llm) ! mass in each cell ! not used
107	REAL, INTENT(IN) :: pq(iip1, jjp1, llm, nqtot) ! tracers
108	REAL, INTENT(IN) :: pphis(iip1, jjp1) ! surface geopotential
109	REAL, INTENT(IN) :: pphi(iip1, jjp1, llm) ! geopotential
110
111	REAL, INTENT(IN) :: pdvcov(iip1, jjm, llm) ! dynamical tendency on vcov
112	REAL, INTENT(IN) :: pducov(iip1, jjp1, llm) ! dynamical tendency on ucov
113	REAL, INTENT(IN) :: pdteta(iip1, jjp1, llm) ! dynamical tendency on teta
114	! NB: pdteta is used only to compute pcvgt which is in fact not used...
115	REAL, INTENT(IN) :: pdq(iip1, jjp1, llm, nqtot) ! dynamical tendency on tracers
116	! NB: pdq is only used to compute pcvgq which is in fact not used...
117
118	REAL, INTENT(IN) :: pps(iip1, jjp1) ! surface pressure (Pa)
119	REAL, INTENT(IN) :: pp(iip1, jjp1, llmp1) ! pressure at mesh interfaces (Pa)
120	REAL, INTENT(IN) :: ppk(iip1, jjp1, llm) ! Exner at mid-layer
121	REAL, INTENT(IN) :: flxw(iip1, jjp1, llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
122
123	! tendencies (in */s) from the physics
124	REAL, INTENT(OUT) :: pdvfi(iip1, jjm, llm) ! tendency on covariant meridional wind
125	REAL, INTENT(OUT) :: pdufi(iip1, jjp1, llm) ! tendency on covariant zonal wind
126	REAL, INTENT(OUT) :: pdhfi(iip1, jjp1, llm) ! tendency on potential temperature (K/s)
127	REAL, INTENT(OUT) :: pdqfi(iip1, jjp1, llm, nqtot) ! tendency on tracers
128	REAL, INTENT(OUT) :: pdpsfi(iip1, jjp1) ! tendency on surface pressure (Pa/s)
129
130
131	! Local variables :
132	! -----------------
133
134	INTEGER :: i, j, l, ig0, ig, iq, itr
135	REAL :: zpsrf(ngridmx)
136	REAL :: zplev(ngridmx, llm + 1), zplay(ngridmx, llm)
137	REAL :: zphi(ngridmx, llm), zphis(ngridmx)
138
139	REAL :: zrot(iip1, jjm, llm) ! AdlC May 2014
140	REAL :: zufi(ngridmx, llm), zvfi(ngridmx, llm)
141	REAL :: zrfi(ngridmx, llm) ! relative wind vorticity
142	REAL :: ztfi(ngridmx, llm), zqfi(ngridmx, llm, nqtot)
143	REAL :: zpk(ngridmx, llm)
144
145	REAL :: pcvgu(ngridmx, llm), pcvgv(ngridmx, llm)
146	REAL :: pcvgt(ngridmx, llm), pcvgq(ngridmx, llm, 2)
147
148	REAL :: zdufi(ngridmx, llm), zdvfi(ngridmx, llm)
149	REAL :: zdtfi(ngridmx, llm), zdqfi(ngridmx, llm, nqtot)
150	REAL :: zdpsrf(ngridmx)
151
152	REAL :: zdufic(ngridmx, llm), zdvfic(ngridmx, llm)
153	REAL :: zdtfic(ngridmx, llm), zdqfic(ngridmx, llm, nqtot)
154	REAL :: jH_cur_split, zdt_split
155	LOGICAL :: debut_split, lafin_split
156	INTEGER :: isplit
157
158	REAL :: zsin(iim), zcos(iim), z1(iim)
159	REAL :: zsinbis(iim), zcosbis(iim), z1bis(iim)
160	REAL :: unskap, pksurcp
161
162	REAL :: flxwfi(ngridmx, llm) ! Flux de masse verticale sur la grille physiq
163
164	LOGICAL, SAVE :: firstcal = .TRUE., debut = .TRUE.
165	! REAL rdayvrai
166
167
168	!-----------------------------------------------------------------------
169
170	! 1. Initialisations :
171	! --------------------
172
173
174	IF (firstcal) THEN
175	debut = .TRUE.
176	IF (ngridmx/=2 + (jjm - 1) * iim) THEN
177	WRITE(lunout, *) 'STOP dans calfis'
178	WRITE(lunout, *) &
179	'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
180	WRITE(lunout, *) ' ngridmx jjm iim '
181	WRITE(lunout, *) ngridmx, jjm, iim
182	CALL abort_gcm("calfis", "", 1)
183	ENDIF
184	ELSE
185	debut = .FALSE.
186	ENDIF ! of IF (firstcal)
187
188
189
190	!-----------------------------------------------------------------------
191	! 40. transformation des variables dynamiques en variables physiques:
192	! ---------------------------------------------------------------
193
194	! 41. pressions au sol (en Pascals)
195	! ----------------------------------
196
197	zpsrf(1) = pps(1, 1)
198
199	ig0 = 2
200	DO j = 2, jjm
201	CALL SCOPY(iim, pps(1, j), 1, zpsrf(ig0), 1)
202	ig0 = ig0 + iim
203	ENDDO
204
205	zpsrf(ngridmx) = pps(1, jjp1)
206
207
208	! 42. pression intercouches et fonction d'Exner:
209
210	! -----------------------------------------------------------------
211	! .... zplev definis aux (llm +1) interfaces des couches ....
212	! .... zplay definis aux ( llm ) milieux des couches ....
213	! -----------------------------------------------------------------
214
215	! ... Exner = cp * ( p(l) / preff ) ** kappa ....
216
217	unskap = 1. / kappa
218
219	DO l = 1, llm
220	zpk(1, l) = ppk(1, 1, l)
221	zplev(1, l) = pp(1, 1, l)
222	ig0 = 2
223	DO j = 2, jjm
224	DO i = 1, iim
225	zpk(ig0, l) = ppk(i, j, l)
226	zplev(ig0, l) = pp(i, j, l)
227	ig0 = ig0 + 1
228	ENDDO
229	ENDDO
230	zpk(ngridmx, l) = ppk(1, jjp1, l)
231	zplev(ngridmx, l) = pp(1, jjp1, l)
232	ENDDO
233	zplev(1, llmp1) = pp(1, 1, llmp1)
234	ig0 = 2
235	DO j = 2, jjm
236	DO i = 1, iim
237	zplev(ig0, llmp1) = pp(i, j, llmp1)
238	ig0 = ig0 + 1
239	ENDDO
240	ENDDO
241	zplev(ngridmx, llmp1) = pp(1, jjp1, llmp1)
242
243	!
244
245	! 43. temperature naturelle (en K) et pressions milieux couches .
246	! ---------------------------------------------------------------
247
248	DO l = 1, llm
249
250	pksurcp = ppk(1, 1, l) / cpp
251	zplay(1, l) = preff * pksurcp ** unskap
252	ztfi(1, l) = pteta(1, 1, l) * pksurcp
253	pcvgt(1, l) = pdteta(1, 1, l) * pksurcp / pmasse(1, 1, l)
254	ig0 = 2
255
256	DO j = 2, jjm
257	DO i = 1, iim
258	pksurcp = ppk(i, j, l) / cpp
259	zplay(ig0, l) = preff * pksurcp ** unskap
260	ztfi(ig0, l) = pteta(i, j, l) * pksurcp
261	pcvgt(ig0, l) = pdteta(i, j, l) * pksurcp / pmasse(i, j, l)
262	ig0 = ig0 + 1
263	ENDDO
264	ENDDO
265
266	pksurcp = ppk(1, jjp1, l) / cpp
267	zplay(ig0, l) = preff * pksurcp ** unskap
268	ztfi (ig0, l) = pteta(1, jjp1, l) * pksurcp
269	pcvgt(ig0, l) = pdteta(1, jjp1, l) * pksurcp / pmasse(1, jjp1, l)
270
271	ENDDO
272
273	! 43.bis traceurs
274	! ---------------
275
276	itr = 0
277	DO iq = 1, nqtot
278	IF(.NOT.tracers(iq)%isAdvected) CYCLE
279	itr = itr + 1
280	DO l = 1, llm
281	zqfi(1, l, itr) = pq(1, 1, l, iq)
282	ig0 = 2
283	DO j = 2, jjm
284	DO i = 1, iim
285	zqfi(ig0, l, itr) = pq(i, j, l, iq)
286	ig0 = ig0 + 1
287	ENDDO
288	ENDDO
289	zqfi(ig0, l, itr) = pq(1, jjp1, l, iq)
290	ENDDO
291	ENDDO
292
293	! convergence dynamique pour les traceurs "EAU"
294	! Earth-specific treatment of first 2 tracers (water)
295	IF (planet_type=="earth") THEN
296	DO iq = 1, 2
297	DO l = 1, llm
298	pcvgq(1, l, iq) = pdq(1, 1, l, iq) / pmasse(1, 1, l)
299	ig0 = 2
300	DO j = 2, jjm
301	DO i = 1, iim
302	pcvgq(ig0, l, iq) = pdq(i, j, l, iq) / pmasse(i, j, l)
303	ig0 = ig0 + 1
304	ENDDO
305	ENDDO
306	pcvgq(ig0, l, iq) = pdq(1, jjp1, l, iq) / pmasse(1, jjp1, l)
307	ENDDO
308	ENDDO
309	endif ! of if (planet_type=="earth")
310
311
312	! Geopotentiel calcule par rapport a la surface locale:
313	! -----------------------------------------------------
314
315	CALL gr_dyn_fi(llm, iip1, jjp1, ngridmx, pphi, zphi)
316	CALL gr_dyn_fi(1, iip1, jjp1, ngridmx, pphis, zphis)
317	DO l = 1, llm
318	DO ig = 1, ngridmx
319	zphi(ig, l) = zphi(ig, l) - zphis(ig)
320	ENDDO
321	ENDDO
322
323	! .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
324	! JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
325	! de masse est calclue dans advtrac.F
326	! DO l=1,llm
327	! pvervel(1,l)=pw(1,1,l) * g /apoln
328	! ig0=2
329	! DO j=2,jjm
330	! DO i = 1, iim
331	! pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
332	! ig0 = ig0 + 1
333	! ENDDO
334	! ENDDO
335	! pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
336	! ENDDO
337
338
339	! 45. champ u:
340	! ------------
341
342	DO l = 1, llm
343
344	DO j = 2, jjm
345	ig0 = 1 + (j - 2) * iim
346	zufi(ig0 + 1, l) = 0.5 * &
347	(pucov(iim, j, l) / cu(iim, j) + pucov(1, j, l) / cu(1, j))
348	pcvgu(ig0 + 1, l) = 0.5 * &
349	(pducov(iim, j, l) / cu(iim, j) + pducov(1, j, l) / cu(1, j))
350	DO i = 2, iim
351	zufi(ig0 + i, l) = 0.5 * &
352	(pucov(i - 1, j, l) / cu(i - 1, j) + pucov(i, j, l) / cu(i, j))
353	pcvgu(ig0 + i, l) = 0.5 * &
354	(pducov(i - 1, j, l) / cu(i - 1, j) + pducov(i, j, l) / cu(i, j))
355	END DO
356	END DO
357
358	END DO
359
360
361	! Alvaro de la Camara (May 2014)
362	! 46.1 Calcul de la vorticite et passage sur la grille physique
363	! --------------------------------------------------------------
364	DO l = 1, llm
365	DO i = 1, iim
366	DO j = 1, jjm
367	zrot(i, j, l) = (pvcov(i + 1, j, l) - pvcov(i, j, l) &
368	+ pucov(i, j + 1, l) - pucov(i, j, l)) &
369	/ (cu(i, j) + cu(i, j + 1)) &
370	/ (cv(i + 1, j) + cv(i, j)) * 4
371	enddo
372	enddo
373	ENDDO
374
375	! 46.champ v:
376	! -----------
377
378	DO l = 1, llm
379	DO j = 2, jjm
380	ig0 = 1 + (j - 2) * iim
381	DO i = 1, iim
382	zvfi(ig0 + i, l) = 0.5 * &
383	(pvcov(i, j - 1, l) / cv(i, j - 1) + pvcov(i, j, l) / cv(i, j))
384	pcvgv(ig0 + i, l) = 0.5 * &
385	(pdvcov(i, j - 1, l) / cv(i, j - 1) + pdvcov(i, j, l) / cv(i, j))
386	ENDDO
387	zrfi(ig0 + 1, l) = 0.25 * (zrot(iim, j - 1, l) + zrot(iim, j, l) &
388	+ zrot(1, j - 1, l) + zrot(1, j, l))
389	DO i = 2, iim
390	zrfi(ig0 + i, l) = 0.25 * (zrot(i - 1, j - 1, l) + zrot(i - 1, j, l) &
391	+ zrot(i, j - 1, l) + zrot(i, j, l)) ! AdlC MAY 2014
392	ENDDO
393	ENDDO
394	ENDDO
395
396
397	! 47. champs de vents aux pole nord
398	! ------------------------------
399	! U = 1 / pi * integrale [ v * cos(long) * d long ]
400	! V = 1 / pi * integrale [ v * sin(long) * d long ]
401
402	DO l = 1, llm
403
404	z1(1) = (rlonu(1) - rlonu(iim) + 2. * pi) * pvcov(1, 1, l) / cv(1, 1)
405	z1bis(1) = (rlonu(1) - rlonu(iim) + 2. * pi) * pdvcov(1, 1, l) / cv(1, 1)
406	DO i = 2, iim
407	z1(i) = (rlonu(i) - rlonu(i - 1)) * pvcov(i, 1, l) / cv(i, 1)
408	z1bis(i) = (rlonu(i) - rlonu(i - 1)) * pdvcov(i, 1, l) / cv(i, 1)
409	ENDDO
410
411	DO i = 1, iim
412	zcos(i) = COS(rlonv(i)) * z1(i)
413	zcosbis(i) = COS(rlonv(i)) * z1bis(i)
414	zsin(i) = SIN(rlonv(i)) * z1(i)
415	zsinbis(i) = SIN(rlonv(i)) * z1bis(i)
416	ENDDO
417
418	zufi(1, l) = SSUM(iim, zcos, 1) / pi
419	pcvgu(1, l) = SSUM(iim, zcosbis, 1) / pi
420	zvfi(1, l) = SSUM(iim, zsin, 1) / pi
421	pcvgv(1, l) = SSUM(iim, zsinbis, 1) / pi
422	zrfi(1, l) = 0.
423	ENDDO
424
425
426	! 48. champs de vents aux pole sud:
427	! ---------------------------------
428	! U = 1 / pi * integrale [ v * cos(long) * d long ]
429	! V = 1 / pi * integrale [ v * sin(long) * d long ]
430
431	DO l = 1, llm
432
433	z1(1) = (rlonu(1) - rlonu(iim) + 2. * pi) * pvcov(1, jjm, l) / cv(1, jjm)
434	z1bis(1) = (rlonu(1) - rlonu(iim) + 2. * pi) * pdvcov(1, jjm, l) / cv(1, jjm)
435	DO i = 2, iim
436	z1(i) = (rlonu(i) - rlonu(i - 1)) * pvcov(i, jjm, l) / cv(i, jjm)
437	z1bis(i) = (rlonu(i) - rlonu(i - 1)) * pdvcov(i, jjm, l) / cv(i, jjm)
438	ENDDO
439
440	DO i = 1, iim
441	zcos(i) = COS(rlonv(i)) * z1(i)
442	zcosbis(i) = COS(rlonv(i)) * z1bis(i)
443	zsin(i) = SIN(rlonv(i)) * z1(i)
444	zsinbis(i) = SIN(rlonv(i)) * z1bis(i)
445	ENDDO
446
447	zufi(ngridmx, l) = SSUM(iim, zcos, 1) / pi
448	pcvgu(ngridmx, l) = SSUM(iim, zcosbis, 1) / pi
449	zvfi(ngridmx, l) = SSUM(iim, zsin, 1) / pi
450	pcvgv(ngridmx, l) = SSUM(iim, zsinbis, 1) / pi
451	zrfi(ngridmx, l) = 0.
452	ENDDO
453
454	! On change de grille, dynamique vers physiq, pour le flux de masse verticale
455	CALL gr_dyn_fi(llm, iip1, jjp1, ngridmx, flxw, flxwfi)
456
457	!-----------------------------------------------------------------------
458	! Appel de la physique:
459	! ---------------------
460
461
462
463	! WRITE(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
464	zdt_split = dtphys / nsplit_phys
465	zdufic(:, :) = 0.
466	zdvfic(:, :) = 0.
467	zdtfic(:, :) = 0.
468	zdqfic(:, :, :) = 0.
469
470	IF (CPPKEY_PHYS) THEN
471
472	DO isplit = 1, nsplit_phys
473
474	jH_cur_split = jH_cur + (isplit - 1) * dtvr / (daysec * nsplit_phys)
475	debut_split = debut.AND.isplit==1
476	lafin_split = lafin.AND.isplit==nsplit_phys
477
478	! if (planet_type=="earth") THEN
479	CALL call_physiq(ngridmx, llm, nqtot, tracers(:)%name, &
480	debut_split, lafin_split, &
481	jD_cur, jH_cur_split, zdt_split, &
482	zplev, zplay, &
483	zpk, zphi, zphis, &
484	presnivs, &
485	zufi, zvfi, zrfi, ztfi, zqfi, &
486	flxwfi, pducov, &
487	zdufi, zdvfi, zdtfi, zdqfi, zdpsrf)
488
489	! ELSE IF ( planet_type=="generic" ) THEN
490	! CALL physiq (ngridmx, !! ngrid
491	! . llm, !! nlayer
492	! . nqtot, !! nq
493	! . tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
494	! . debut_split, !! firstcall
495	! . lafin_split, !! lastcall
496	! . jD_cur, !! pday. see leapfrog
497	! . jH_cur_split, !! ptime "fraction of day"
498	! . zdt_split, !! ptimestep
499	! . zplev, !! pplev
500	! . zplay, !! pplay
501	! . zphi, !! pphi
502	! . zufi, !! pu
503	! . zvfi, !! pv
504	! . ztfi, !! pt
505	! . zqfi, !! pq
506	! . flxwfi, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
507	! . zdufi, !! pdu
508	! . zdvfi, !! pdv
509	! . zdtfi, !! pdt
510	! . zdqfi, !! pdq
511	! . zdpsrf, !! pdpsrf
512	! . tracerdyn) !! tracerdyn <-- utilite ???
513
514	! ENDIF ! of if (planet_type=="earth")
515
516	zufi(:, :) = zufi(:, :) + zdufi(:, :) * zdt_split
517	zvfi(:, :) = zvfi(:, :) + zdvfi(:, :) * zdt_split
518	ztfi(:, :) = ztfi(:, :) + zdtfi(:, :) * zdt_split
519	zqfi(:, :, :) = zqfi(:, :, :) + zdqfi(:, :, :) * zdt_split
520
521	zdufic(:, :) = zdufic(:, :) + zdufi(:, :)
522	zdvfic(:, :) = zdvfic(:, :) + zdvfi(:, :)
523	zdtfic(:, :) = zdtfic(:, :) + zdtfi(:, :)
524	zdqfic(:, :, :) = zdqfic(:, :, :) + zdqfi(:, :, :)
525
526	enddo ! of do isplit=1,nsplit_phys
527
528	END IF
529
530	zdufi(:, :) = zdufic(:, :) / nsplit_phys
531	zdvfi(:, :) = zdvfic(:, :) / nsplit_phys
532	zdtfi(:, :) = zdtfic(:, :) / nsplit_phys
533	zdqfi(:, :, :) = zdqfic(:, :, :) / nsplit_phys
534
535	!-----------------------------------------------------------------------
536	! transformation des tendances physiques en tendances dynamiques:
537	! ---------------------------------------------------------------
538
539	! tendance sur la pression :
540	! -----------------------------------
541
542	CALL gr_fi_dyn(1, ngridmx, iip1, jjp1, zdpsrf, pdpsfi)
543
544	! 62. enthalpie potentielle
545	! ---------------------
546
547	DO l = 1, llm
548
549	DO i = 1, iip1
550	pdhfi(i, 1, l) = cpp * zdtfi(1, l) / ppk(i, 1, l)
551	pdhfi(i, jjp1, l) = cpp * zdtfi(ngridmx, l) / ppk(i, jjp1, l)
552	ENDDO
553
554	DO j = 2, jjm
555	ig0 = 1 + (j - 2) * iim
556	DO i = 1, iim
557	pdhfi(i, j, l) = cpp * zdtfi(ig0 + i, l) / ppk(i, j, l)
558	ENDDO
559	pdhfi(iip1, j, l) = pdhfi(1, j, l)
560	ENDDO
561
562	ENDDO
563
564
565	! 62. humidite specifique
566	! ---------------------
567	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
568	! DO iq=1,nqtot
569	! DO l=1,llm
570	! DO i=1,iip1
571	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
572	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
573	! ENDDO
574	! DO j=2,jjm
575	! ig0=1+(j-2)*iim
576	! DO i=1,iim
577	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
578	! ENDDO
579	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
580	! ENDDO
581	! ENDDO
582	! ENDDO
583
584	! 63. traceurs
585	! ------------
586	! initialisation des tendances
587	pdqfi(:, :, :, :) = 0.
588
589	itr = 0
590	DO iq = 1, nqtot
591	IF(.NOT.tracers(iq)%isAdvected) CYCLE
592	itr = itr + 1
593	DO l = 1, llm
594	DO i = 1, iip1
595	pdqfi(i, 1, l, iq) = zdqfi(1, l, itr)
596	pdqfi(i, jjp1, l, iq) = zdqfi(ngridmx, l, itr)
597	ENDDO
598	DO j = 2, jjm
599	ig0 = 1 + (j - 2) * iim
600	DO i = 1, iim
601	pdqfi(i, j, l, iq) = zdqfi(ig0 + i, l, itr)
602	ENDDO
603	pdqfi(iip1, j, l, iq) = pdqfi(1, j, l, itr)
604	ENDDO
605	ENDDO
606	ENDDO
607
608	! 65. champ u:
609	! ------------
610
611	DO l = 1, llm
612
613	DO i = 1, iip1
614	pdufi(i, 1, l) = 0.
615	pdufi(i, jjp1, l) = 0.
616	ENDDO
617
618	DO j = 2, jjm
619	ig0 = 1 + (j - 2) * iim
620	DO i = 1, iim - 1
621	pdufi(i, j, l) = &
622	0.5 * (zdufi(ig0 + i, l) + zdufi(ig0 + i + 1, l)) * cu(i, j)
623	ENDDO
624	pdufi(iim, j, l) = &
625	0.5 * (zdufi(ig0 + 1, l) + zdufi(ig0 + iim, l)) * cu(iim, j)
626	pdufi(iip1, j, l) = pdufi(1, j, l)
627	ENDDO
628
629	ENDDO
630
631
632	! 67. champ v:
633	! ------------
634
635	DO l = 1, llm
636
637	DO j = 2, jjm - 1
638	ig0 = 1 + (j - 2) * iim
639	DO i = 1, iim
640	pdvfi(i, j, l) = &
641	0.5 * (zdvfi(ig0 + i, l) + zdvfi(ig0 + i + iim, l)) * cv(i, j)
642	ENDDO
643	pdvfi(iip1, j, l) = pdvfi(1, j, l)
644	ENDDO
645	ENDDO
646
647
648	! 68. champ v pres des poles:
649	! ---------------------------
650	! v = U * cos(long) + V * SIN(long)
651
652	DO l = 1, llm
653
654	DO i = 1, iim
655	pdvfi(i, 1, l) = &
656	zdufi(1, l) * COS(rlonv(i)) + zdvfi(1, l) * SIN(rlonv(i))
657	pdvfi(i, jjm, l) = zdufi(ngridmx, l) * COS(rlonv(i)) &
658	+ zdvfi(ngridmx, l) * SIN(rlonv(i))
659	pdvfi(i, 1, l) = &
660	0.5 * (pdvfi(i, 1, l) + zdvfi(i + 1, l)) * cv(i, 1)
661	pdvfi(i, jjm, l) = &
662	0.5 * (pdvfi(i, jjm, l) + zdvfi(ngridmx - iip1 + i, l)) * cv(i, jjm)
663	ENDDO
664
665	pdvfi(iip1, 1, l) = pdvfi(1, 1, l)
666	pdvfi(iip1, jjm, l) = pdvfi(1, jjm, l)
667
668	ENDDO
669
670	!-----------------------------------------------------------------------
671	firstcal = .FALSE.
672
673	RETURN
674	END SUBROUTINE calfis

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

Download in other formats:

Original Format