Context Navigation

calfis.F @ 4283

Last change on this file since 4283 was 4056, checked in by dcugnet, 2 years ago

Most of the changes are intended to help to eventually remove the constraints about the tracers assumptions, in particular water tracers.

Remove index tables itr_indice and niadv, replaced by tracers(:)%isAdvected and tracers(:)%isH2OFamily. Most of the loops are now from 1 to nqtot:
- DO iq=nqo+1,nqtot loops are replaced with: DO iq=1,nqtot

IF(tracers(iq)%isH2Ofamily) CYCLE

DO it=1,nbtr; iq=niadv(it+nqo)

and DO it=1,nqtottr; iq=itr_indice(it) loops are replaced with:

it = 0
DO iq = 1, nqtot

IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
it = it+1

Move some StratAer? related code from infotrac to infotrac_phy
Remove "nqperes" variable:

DO iq=1,nqpere loops are replaced with:
DO iq=1,nqtot

IF(tracers(iq)%parent/='air') CYCLE

Cosmetic changes (justification, SELECT CASE instead of multiple IF...) mostly in advtrac* routines.

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: 21.2 KB

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

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

Download in other formats:

Original Format