Context Navigation

calfis.f90 @ 5174

Last change on this file since 5174 was 5159, checked in by abarral, 3 months ago
Put dimensions.h and paramet.h into 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 5159 2024-08-02 19:58:25Z 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 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