Context Navigation

calfis.f90 @ 5133

Last change on this file since 5133 was 5128, checked in by abarral, 5 months ago
Correct bug in vlspltqs_loc.f90 from r2270 where we call SSUM with incorrect arguments. Merge the three different versions of abort_gcm into one Fix seq, para 3D compilation broken from r5107 onwards (lint) usual + Remove uneeded fixed-form continuations
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.4 KB

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

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

Download in other formats:

Original Format