source: LMDZ6/trunk/libf/dynphy_lonlat/calfis.F90 @ 5254

Last change on this file since 5254 was 5250, checked in by abarral, 6 hours ago

Wrap uses of cpp key CPP_PHYS

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