Context Navigation

calfis.f90 @ 5137

Last change on this file since 5137 was 5136, checked in by abarral, 3 months ago
Put comgeom.h, comgeom2.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.4 KB

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

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

Download in other formats:

Original Format