source: LMDZ6/trunk/libf/dynphy_lonlat/calfis.f90 @ 5272

Last change on this file since 5272 was 5272, checked in by abarral, 8 weeks ago

Turn paramet.h into a module

  • 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.9 KB
Line 
1!
2! $Id: calfis.f90 5272 2024-10-24 15:53:15Z 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  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
38  USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
39          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
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  INTEGER :: ngridmx
94  PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
95
96  include "comgeom2.h"
97  include "iniprint.h"
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
165  REAL :: SSUM
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.NE.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
201  zpsrf(1) = pps(1,1)
202
203  ig0  = 2
204  DO j = 2,jjm
205     CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
206     ig0 = ig0+iim
207  ENDDO
208
209  zpsrf(ngridmx) = pps(1,jjp1)
210
211
212  !   42. pression intercouches et fonction d'Exner:
213  !
214  !   -----------------------------------------------------------------
215  ! .... zplev  definis aux (llm +1) interfaces des couches  ....
216  ! .... zplay  definis aux (  llm )    milieux des couches  ....
217  !   -----------------------------------------------------------------
218
219  !    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
220  !
221   unskap   = 1./ kappa
222  !
223  DO l = 1, llm
224    zpk(   1,l ) = ppk(1,1,l)
225    zplev( 1,l ) = pp(1,1,l)
226    ig0 = 2
227      DO j = 2, jjm
228         DO i =1, iim
229          zpk(   ig0,l ) = ppk(i,j,l)
230          zplev( ig0,l ) = pp(i,j,l)
231          ig0 = ig0 +1
232         ENDDO
233      ENDDO
234    zpk(   ngridmx,l ) = ppk(1,jjp1,l)
235    zplev( ngridmx,l ) = pp(1,jjp1,l)
236  ENDDO
237    zplev( 1,llmp1 ) = pp(1,1,llmp1)
238    ig0 = 2
239      DO j = 2, jjm
240         DO i =1, iim
241          zplev( ig0,llmp1 ) = pp(i,j,llmp1)
242          ig0 = ig0 +1
243         ENDDO
244      ENDDO
245    zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
246  !
247  !
248
249  !   43. temperature naturelle (en K) et pressions milieux couches .
250  !   ---------------------------------------------------------------
251
252  DO l=1,llm
253
254     pksurcp     =  ppk(1,1,l) / cpp
255     zplay(1,l)  =  preff * pksurcp ** unskap
256     ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
257     pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
258     ig0         = 2
259
260     DO j = 2, jjm
261        DO i = 1, iim
262          pksurcp        = ppk(i,j,l) / cpp
263          zplay(ig0,l)   = preff * pksurcp ** unskap
264          ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
265          pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
266          ig0            = ig0 + 1
267        ENDDO
268     ENDDO
269
270     pksurcp       = ppk(1,jjp1,l) / cpp
271     zplay(ig0,l)  = preff * pksurcp ** unskap
272     ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
273     pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
274
275  ENDDO
276
277  !   43.bis traceurs
278  !   ---------------
279  !
280  itr=0
281  DO iq=1,nqtot
282     IF(.NOT.tracers(iq)%isAdvected) CYCLE
283     itr = itr + 1
284     DO l=1,llm
285        zqfi(1,l,itr) = pq(1,1,l,iq)
286        ig0           = 2
287        DO j=2,jjm
288           DO i = 1, iim
289              zqfi(ig0,l,itr) = pq(i,j,l,iq)
290              ig0             = ig0 + 1
291           ENDDO
292        ENDDO
293        zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
294     ENDDO
295  ENDDO
296
297  !   convergence dynamique pour les traceurs "EAU"
298  ! Earth-specific treatment of first 2 tracers (water)
299   if (planet_type=="earth") then
300    DO iq=1,2
301     DO l=1,llm
302        pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
303        ig0          = 2
304        DO j=2,jjm
305           DO i = 1, iim
306              pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
307              ig0             = ig0 + 1
308           ENDDO
309        ENDDO
310        pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
311     ENDDO
312    ENDDO
313   endif ! of if (planet_type=="earth")
314
315
316  !   Geopotentiel calcule par rapport a la surface locale:
317  !   -----------------------------------------------------
318
319  CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
320  CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
321  DO l=1,llm
322     DO ig=1,ngridmx
323       zphi(ig,l)=zphi(ig,l)-zphis(ig)
324     ENDDO
325  ENDDO
326
327  !   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
328  ! JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
329  !    de masse est calclue dans advtrac.F
330   ! DO l=1,llm
331   !   pvervel(1,l)=pw(1,1,l) * g /apoln
332   !   ig0=2
333   !  DO j=2,jjm
334   !      DO i = 1, iim
335   !         pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
336   !         ig0 = ig0 + 1
337   !      ENDDO
338  !       ENDDO
339   !   pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
340   ! ENDDO
341
342  !
343  !   45. champ u:
344  !   ------------
345
346  DO l=1,llm
347
348     DO j=2,jjm
349        ig0 = 1+(j-2)*iim
350        zufi(ig0+1,l)= 0.5 * &
351              ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
352        pcvgu(ig0+1,l)= 0.5 * &
353              ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
354        DO i=2,iim
355           zufi(ig0+i,l)= 0.5 * &
356                 ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
357           pcvgu(ig0+i,l)= 0.5 * &
358                 ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
359        END DO
360     END DO
361
362  END DO
363
364
365  !  Alvaro de la Camara (May 2014)
366  !  46.1 Calcul de la vorticite et passage sur la grille physique
367  !  --------------------------------------------------------------
368  DO l=1,llm
369    do i=1,iim
370      do j=1,jjm
371        zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l) &
372              + pucov(i,j+1,l) - pucov(i,j,l)) &
373              / (cu(i,j)+cu(i,j+1)) &
374              / (cv(i+1,j)+cv(i,j)) *4
375      enddo
376    enddo
377  ENDDO
378
379  !   46.champ v:
380  !   -----------
381
382  DO l=1,llm
383     DO j=2,jjm
384        ig0=1+(j-2)*iim
385        DO i=1,iim
386           zvfi(ig0+i,l)= 0.5 * &
387                 ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
388           pcvgv(ig0+i,l)= 0.5 * &
389                 ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
390        ENDDO
391           zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l) &
392                 +zrot(1,j-1,l)+zrot(1,j,l))
393        DO i=2,iim
394           zrfi(ig0 + i,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l) &
395                 +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
396        ENDDO
397     ENDDO
398  ENDDO
399
400
401  !   47. champs de vents aux pole nord
402  !   ------------------------------
403     ! U = 1 / pi  *  integrale [ v * cos(long) * d long ]
404     ! V = 1 / pi  *  integrale [ v * sin(long) * d long ]
405
406  DO l=1,llm
407
408     z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
409     z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
410     DO i=2,iim
411        z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
412        z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
413     ENDDO
414
415     DO i=1,iim
416        zcos(i)   = COS(rlonv(i))*z1(i)
417        zcosbis(i)= COS(rlonv(i))*z1bis(i)
418        zsin(i)   = SIN(rlonv(i))*z1(i)
419        zsinbis(i)= SIN(rlonv(i))*z1bis(i)
420     ENDDO
421
422     zufi(1,l)  = SSUM(iim,zcos,1)/pi
423     pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
424     zvfi(1,l)  = SSUM(iim,zsin,1)/pi
425     pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
426     zrfi(1, l) = 0.
427  ENDDO
428
429
430  !   48. champs de vents aux pole sud:
431  !   ---------------------------------
432     ! U = 1 / pi  *  integrale [ v * cos(long) * d long ]
433     ! V = 1 / pi  *  integrale [ v * sin(long) * d long ]
434
435  DO l=1,llm
436
437     z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
438     z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
439     DO i=2,iim
440        z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
441        z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
442     ENDDO
443
444     DO i=1,iim
445        zcos(i)    = COS(rlonv(i))*z1(i)
446        zcosbis(i) = COS(rlonv(i))*z1bis(i)
447        zsin(i)    = SIN(rlonv(i))*z1(i)
448        zsinbis(i) = SIN(rlonv(i))*z1bis(i)
449     ENDDO
450
451     zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
452     pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
453     zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
454     pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
455     zrfi(ngridmx, l) = 0.
456  ENDDO
457  !
458  ! On change de grille, dynamique vers physiq, pour le flux de masse verticale
459  CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
460
461  !-----------------------------------------------------------------------
462  !   Appel de la physique:
463  !   ---------------------
464
465
466
467   ! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
468  zdt_split=dtphys/nsplit_phys
469  zdufic(:,:)=0.
470  zdvfic(:,:)=0.
471  zdtfic(:,:)=0.
472  zdqfic(:,:,:)=0.
473
474IF (CPPKEY_PHYS) THEN
475
476   do isplit=1,nsplit_phys
477
478     jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
479     debut_split=debut.and.isplit==1
480     lafin_split=lafin.and.isplit==nsplit_phys
481
482   ! if (planet_type=="earth") then
483    CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name, &
484          debut_split,lafin_split, &
485          jD_cur,jH_cur_split,zdt_split, &
486          zplev,zplay, &
487          zpk,zphi,zphis, &
488          presnivs, &
489          zufi,zvfi,zrfi,ztfi,zqfi, &
490          flxwfi,pducov, &
491          zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
492  !
493  !  else if ( planet_type=="generic" ) then
494  !
495  !     CALL physiq (ngridmx,     !! ngrid
496  ! .             llm,            !! nlayer
497  ! .             nqtot,          !! nq
498  ! .             tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
499  ! .             debut_split,    !! firstcall
500  ! .             lafin_split,    !! lastcall
501  ! .             jD_cur,         !! pday. see leapfrog
502  ! .             jH_cur_split,   !! ptime "fraction of day"
503  ! .             zdt_split,      !! ptimestep
504  ! .             zplev,          !! pplev
505  ! .             zplay,          !! pplay
506  ! .             zphi,           !! pphi
507  ! .             zufi,           !! pu
508  ! .             zvfi,           !! pv
509  ! .             ztfi,           !! pt
510  ! .             zqfi,           !! pq
511  ! .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
512  ! .             zdufi,          !! pdu
513  ! .             zdvfi,          !! pdv
514  ! .             zdtfi,          !! pdt
515  ! .             zdqfi,          !! pdq
516  ! .             zdpsrf,         !! pdpsrf
517  ! .             tracerdyn)      !! tracerdyn <-- utilite ???
518  !
519  !  endif ! of if (planet_type=="earth")
520
521     zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
522     zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
523     ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
524     zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
525
526     zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
527     zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
528     zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
529     zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
530
531   enddo ! of do isplit=1,nsplit_phys
532
533END IF
534  ! of #ifdef CPP_PHYS
535
536  zdufi(:,:)=zdufic(:,:)/nsplit_phys
537  zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
538  zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
539  zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
540
541
542500   CONTINUE
543
544  !-----------------------------------------------------------------------
545  !   transformation des tendances physiques en tendances dynamiques:
546  !   ---------------------------------------------------------------
547
548  !  tendance sur la pression :
549  !  -----------------------------------
550
551  CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
552  !
553  !   62. enthalpie potentielle
554  !   ---------------------
555
556  DO l=1,llm
557
558     DO i=1,iip1
559      pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
560      pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
561     ENDDO
562
563     DO j=2,jjm
564        ig0=1+(j-2)*iim
565        DO i=1,iim
566           pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
567        ENDDO
568           pdhfi(iip1,j,l) =  pdhfi(1,j,l)
569     ENDDO
570
571  ENDDO
572
573
574  !   62. humidite specifique
575  !   ---------------------
576  ! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
577   ! DO iq=1,nqtot
578   !    DO l=1,llm
579   !       DO i=1,iip1
580   !          pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
581   !          pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
582   !       ENDDO
583   !       DO j=2,jjm
584   !          ig0=1+(j-2)*iim
585   !          DO i=1,iim
586   !             pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
587   !          ENDDO
588   !          pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
589   !       ENDDO
590   !    ENDDO
591   ! ENDDO
592
593  !   63. traceurs
594  !   ------------
595  ! initialisation des tendances
596  pdqfi(:,:,:,:)=0.
597  !
598  itr = 0
599  DO iq=1,nqtot
600     IF(.NOT.tracers(iq)%isAdvected) CYCLE
601     itr = itr + 1
602     DO l=1,llm
603        DO i=1,iip1
604           pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
605           pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
606        ENDDO
607        DO j=2,jjm
608           ig0=1+(j-2)*iim
609           DO i=1,iim
610              pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
611           ENDDO
612           pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
613        ENDDO
614     ENDDO
615  ENDDO
616
617  !   65. champ u:
618  !   ------------
619
620  DO l=1,llm
621
622     DO i=1,iip1
623        pdufi(i,1,l)    = 0.
624        pdufi(i,jjp1,l) = 0.
625     ENDDO
626
627     DO j=2,jjm
628        ig0=1+(j-2)*iim
629        DO i=1,iim-1
630           pdufi(i,j,l)= &
631                 0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
632        ENDDO
633        pdufi(iim,j,l)= &
634              0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
635        pdufi(iip1,j,l)=pdufi(1,j,l)
636     ENDDO
637
638  ENDDO
639
640
641  !   67. champ v:
642  !   ------------
643
644  DO l=1,llm
645
646     DO j=2,jjm-1
647        ig0=1+(j-2)*iim
648        DO i=1,iim
649           pdvfi(i,j,l)= &
650                 0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
651        ENDDO
652        pdvfi(iip1,j,l) = pdvfi(1,j,l)
653     ENDDO
654  ENDDO
655
656
657  !   68. champ v pres des poles:
658  !   ---------------------------
659   ! v = U * cos(long) + V * SIN(long)
660
661  DO l=1,llm
662
663     DO i=1,iim
664        pdvfi(i,1,l)= &
665              zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
666        pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i)) &
667              +zdvfi(ngridmx,l)*SIN(rlonv(i))
668        pdvfi(i,1,l)= &
669              0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
670        pdvfi(i,jjm,l)= &
671              0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
672      ENDDO
673
674     pdvfi(iip1,1,l)  = pdvfi(1,1,l)
675     pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
676
677  ENDDO
678
679  !-----------------------------------------------------------------------
680
681700   CONTINUE
682
683  firstcal = .FALSE.
684
685  RETURN
686END SUBROUTINE calfis
Note: See TracBrowser for help on using the repository browser.