source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3d/calfis.F

Last change on this file was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

File size: 29.4 KB
Line 
1!
2! $Id: calfis.F 1407 2010-07-07 10:31:52Z fairhead $
3!
4C
5C
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)
28c
29c    Auteur :  P. Le Van, F. Hourdin 
30c   .........
31      USE infotrac, ONLY: nqtot, niadv, tname
32      USE control_mod, ONLY: planet_type, nsplit_phys
33      USE write_field
34      USE cpdet_mod, only: t2tpot,tpot2t
35 
36! used only for zonal averages
37      USE moyzon_mod
38
39      IMPLICIT NONE
40c=======================================================================
41c
42c   1. rearrangement des tableaux et transformation
43c      variables dynamiques  >  variables physiques
44c   2. calcul des termes physiques
45c   3. retransformation des tendances physiques en tendances dynamiques
46c
47c   remarques:
48c   ----------
49c
50c    - les vents sont donnes dans la physique par leurs composantes 
51c      naturelles.
52c    - la variable thermodynamique de la physique est une variable
53c      intensive :   T 
54c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
55c    - les deux seules variables dependant de la geometrie necessaires
56c      pour la physique sont la latitude pour le rayonnement et 
57c      l'aire de la maille quand on veut integrer une grandeur
58c      horizontalement.
59c    - les points de la physique sont les points scalaires de la
60c      la dynamique; numerotation:
61c          1 pour le pole nord
62c          (jjm-1)*iim pour l'interieur du domaine
63c          ngridmx pour le pole sud
64c      ---> ngridmx=2+(jjm-1)*iim
65c
66c     Input :
67c     -------
68c       pucov           covariant zonal velocity
69c       pvcov           covariant meridional velocity 
70c       pteta           potential temperature
71c       pps             surface pressure
72c       pmasse          masse d'air dans chaque maille
73c       pts             surface temperature  (K)
74c       callrad         clef d'appel au rayonnement
75c
76c    Output :
77c    --------
78c        pdufi          tendency for the natural zonal velocity (ms-1)
79c        pdvfi          tendency for the natural meridional velocity 
80c        pdhfi          tendency for the potential temperature (K/s)
81c        pdtsfi         tendency for the surface temperature
82c
83c        pdtrad         radiative tendencies  \  both input
84c        pfluxrad       radiative fluxes      /  and output
85c
86c=======================================================================
87c
88c-----------------------------------------------------------------------
89c
90c    0.  Declarations :
91c    ------------------
92
93#include "dimensions.h"
94#include "paramet.h"
95#include "temps.h"
96#include "logic.h"
97
98      INTEGER ngridmx
99      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
100
101#include "comconst.h"
102#include "comvert.h"
103#include "comgeom2.h"
104#include "iniprint.h"
105
106c    Arguments :
107c    -----------
108      LOGICAL,INTENT(IN) ::  lafin ! .true. for the very last call to physics
109      REAL,INTENT(IN) :: jD_cur, jH_cur
110      REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
111      REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
112      REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
113      REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
114      REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
115      REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
116      REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
117
118      REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
119      REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
120      REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
121! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
122! qui lui meme ne sert a rien dans la routine telle qu'elle est
123! ecrite, et que j'ai donc commente....
124      REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
125      ! NB: pdq is only used to compute pcvgq which is in fact not used...
126
127      REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
128      REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
129      REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
130      REAL,INTENT(IN) :: flxw(iip1,jjp1,llm)  ! Vertical mass flux on dynamics grid
131
132      ! tendencies (in */s) from the physics
133      REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
134      REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
135      REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
136      REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
137      REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
138
139
140c    Local variables :
141c    -----------------
142
143      INTEGER i,j,l,ig0,ig,iq,iiq
144      REAL zpsrf(ngridmx)
145      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
146      REAL zphi(ngridmx,llm),zphis(ngridmx)
147
148      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
149      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
150! ADAPTATION GCM POUR CP(T)
151      REAL zteta(ngridmx,llm)
152      REAL zpk(ngridmx,llm)
153
154! RQ SL 13/10/10:
155! Ces calculs ne servent pas.
156! Si necessaire, decommenter ces variables et les calculs...
157!      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
158!      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
159
160      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
161      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
162      REAL zdpsrf(ngridmx)
163
164      REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
165      REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
166      REAL jH_cur_split,zdt_split
167      LOGICAL debut_split,lafin_split
168      INTEGER isplit
169
170      REAL zsin(iim),zcos(iim),z1(iim)
171      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
172      REAL unskap, pksurcp
173      save unskap
174
175      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
176     
177      REAL SSUM
178
179      LOGICAL,SAVE :: firstcal=.true., debut=.true.
180!      REAL rdayvrai
181
182      LOGICAL tracerdyn ! for generic/mars physics call ; possibly to get rid of
183
184! For Titan only right now:
185! to allow for 2D computation of microphys and chemistry
186      LOGICAL,save :: flag_moyzon
187      REAL,allocatable,save :: tmpvar(:,:)
188      REAL,allocatable,save :: tmpvarp1(:,:)
189      REAL,allocatable,save :: tmpvarbar(:)
190      REAL,allocatable,save :: tmpvarbarp1(:)
191      real :: zz1,zz2
192
193c-----------------------------------------------------------------------
194
195c    1. Initialisations :
196c    --------------------
197
198
199      IF ( firstcal )  THEN
200        debut = .TRUE.
201        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
202         write(lunout,*) 'STOP dans calfis'
203         write(lunout,*)
204     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
205         write(lunout,*) '  ngridmx  jjm   iim   '
206         write(lunout,*) ngridmx,jjm,iim
207         STOP
208        ENDIF
209
210        unskap   = 1./ kappa
211
212        flag_moyzon = .false.
213        if(moyzon_ch.or.moyzon_mu) then
214         flag_moyzon = .true.
215         allocate(tmpvar(iip1,llm))
216         allocate(tmpvarp1(iip1,llmp1))
217         allocate(tmpvarbar(llm))
218         allocate(tmpvarbarp1(llmp1))
219        endif
220
221        if (flag_moyzon) call moyzon_init(ngridmx,llm,nqtot)
222
223c----------------------------------------------
224c moyennes globales pour le profil de pression
225       if(planet_type.eq."titan".or.planet_type.eq."venus") then
226        ALLOCATE(plevmoy(llm+1))
227        ALLOCATE(playmoy(llm))
228        ALLOCATE(tmoy(llm))
229        ALLOCATE(tetamoy(llm))
230        ALLOCATE(pkmoy(llm))
231        ALLOCATE(phimoy(0:llm))
232        ALLOCATE(zlevmoy(llm+1))
233        ALLOCATE(zlaymoy(llm))
234        plevmoy=0.
235        do l=1,llmp1
236         do i=1,iip1
237          do j=1,jjp1
238            plevmoy(l)=plevmoy(l)+pp(i,j,l)/(iip1*jjp1)
239          enddo
240         enddo
241        enddo
242        tetamoy=0.
243        pkmoy=0.
244        phimoy=0.
245        do i=1,iip1
246         do j=1,jjp1
247            phimoy(0)=phimoy(0)+pphis(i,j)/(iip1*jjp1)
248         enddo
249        enddo
250        do l=1,llm
251         do i=1,iip1
252          do j=1,jjp1
253            tetamoy(l)=tetamoy(l)+pteta(i,j,l)/(iip1*jjp1)
254            pkmoy(l)=pkmoy(l)+ppk(i,j,l)/(iip1*jjp1)
255            phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1)
256          enddo
257         enddo
258        enddo
259        playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap
260        call tpot2t(llm,tetamoy,tmoy,pkmoy)
261c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
262        zlaymoy(:) = g*rad*rad/(g*rad-phimoy(:))-rad
263        zlevmoy(1) = phimoy(0)/g
264        DO l=2,llm
265            zz1=(playmoy(l-1)+plevmoy(l))/(playmoy(l-1)-plevmoy(l))
266            zz2=(plevmoy(l)  +playmoy(l))/(plevmoy(l)  -playmoy(l))
267            zlevmoy(l)=(zz1*zlaymoy(l-1)+zz2*zlaymoy(l))/(zz1+zz2)
268        ENDDO
269        zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm))
270c-------------------
271c + lat index
272        allocate(klat(ngridmx))
273        klat=0
274        klat(1)  = 1
275        ig0  = 2
276        DO j = 2,jjm
277         do i=0,iim-1
278          klat(ig0+i) = j
279         enddo
280         ig0 = ig0+iim
281        ENDDO
282        klat(ngridmx)  = jjp1
283       endif   ! planet_type=titan
284c----------------------------------------------
285      ELSE
286        debut = .FALSE.
287      ENDIF ! of IF (firstcal)
288
289
290c-----------------------------------------------------------------------
291c   40. transformation des variables dynamiques en variables physiques:
292c   ---------------------------------------------------------------
293
294c   41. pressions au sol (en Pascals)
295c   ----------------------------------
296       
297      zpsrf(1) = pps(1,1)
298
299      ig0  = 2
300      DO j = 2,jjm
301         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
302         ig0 = ig0+iim
303      ENDDO
304
305      zpsrf(ngridmx) = pps(1,jjp1)
306
307c   42. pression intercouches et fonction d'Exner:
308
309c   -----------------------------------------------------------------
310c     .... zplev  definis aux (llm +1) interfaces des couches  ....
311c     .... zplay  definis aux (  llm )    milieux des couches  .... 
312c   -----------------------------------------------------------------
313
314c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
315
316! ADAPTATION GCM POUR CP(T)
317      DO l = 1, llm
318        zpk(   1,l ) = ppk(1,1,l)
319        zteta( 1,l ) = pteta(1,1,l)
320        zplev( 1,l ) = pp(1,1,l)
321        ig0 = 2
322          DO j = 2, jjm
323             DO i =1, iim
324              zpk(   ig0,l ) = ppk(i,j,l)
325              zteta( ig0,l ) = pteta(i,j,l)
326              zplev( ig0,l ) = pp(i,j,l)
327              ig0 = ig0 +1
328             ENDDO
329          ENDDO
330        zpk(   ngridmx,l ) = ppk(1,jjp1,l)
331        zteta( ngridmx,l ) = pteta(1,jjp1,l)
332        zplev( ngridmx,l ) = pp(1,jjp1,l)
333      ENDDO
334        zplev( 1,llmp1 ) = pp(1,1,llmp1)
335        ig0 = 2
336          DO j = 2, jjm
337             DO i =1, iim
338              zplev( ig0,llmp1 ) = pp(i,j,llmp1)
339              ig0 = ig0 +1
340             ENDDO
341          ENDDO
342        zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
343
344      if (flag_moyzon) then
345        tmpvarp1(:,:) = pp(:,1,:)
346        call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
347        zplevbar(1,:) = tmpvarbarp1
348        tmpvar(:,:) = ppk(:,1,:)
349        call moyzon(llm,tmpvar,tmpvarbar)
350        zpkbar(1,:) = tmpvarbar
351        tmpvar(:,:) = pteta(:,1,:)
352        call moyzon(llm,tmpvar,tmpvarbar)
353        ztetabar(1,:) = tmpvarbar
354        call tpot2t(llm,ztetabar(1,:),ztfibar(1,:),zpkbar(1,:))
355        ig0 = 2
356         do j = 2, jjm
357          tmpvarp1(:,:) = pp(:,j,:)
358          call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
359          zplevbar(ig0,:) = tmpvarbarp1
360          tmpvar(:,:) = ppk(:,j,:)
361          call moyzon(llm,tmpvar,tmpvarbar)
362          zpkbar(ig0,:) = tmpvarbar
363          tmpvar(:,:) = pteta(:,j,:)
364          call moyzon(llm,tmpvar,tmpvarbar)
365          ztetabar(ig0,:) = tmpvarbar
366          call tpot2t(llm,ztetabar(ig0,:),ztfibar(ig0,:),zpkbar(ig0,:))
367          ig0 = ig0+1
368          do i=2,iim
369            zplevbar(ig0,:) = zplevbar(ig0-1,:)
370            zpkbar(ig0,:)   = zpkbar(ig0-1,:)
371            ztetabar(ig0,:) = ztetabar(ig0-1,:)
372            ztfibar(ig0,:)  = ztfibar(ig0-1,:)
373            ig0 = ig0+1
374          enddo
375         enddo
376        tmpvarp1(:,:) = pp(:,jjp1,:)
377        call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
378        zplevbar(ngridmx,:) = tmpvarbarp1
379        tmpvar(:,:) = ppk(:,jjp1,:)
380        call moyzon(llm,tmpvar,tmpvarbar)
381        zpkbar(ngridmx,:) = tmpvarbar
382        tmpvar(:,:) = pteta(:,jjp1,:)
383        call moyzon(llm,tmpvar,tmpvarbar)
384        ztetabar(ngridmx,:) = tmpvarbar
385        call tpot2t(llm,ztetabar(ngridmx,:),
386     .                  ztfibar(ngridmx,:),zpkbar(ngridmx,:))
387      endif
388
389c   43. temperature naturelle (en K) et pressions milieux couches .
390c   ---------------------------------------------------------------
391
392! ADAPTATION GCM POUR CP(T)
393         call tpot2t(ngridmx*llm,zteta,ztfi,zpk)
394
395      DO l=1,llm
396
397         pksurcp     =  ppk(1,1,l) / cpp
398         zplay(1,l)  =  preff * pksurcp ** unskap
399!         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
400         ig0         = 2
401
402         DO j = 2, jjm
403            DO i = 1, iim
404              pksurcp        = ppk(i,j,l) / cpp
405              zplay(ig0,l)   = preff * pksurcp ** unskap
406!              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
407              ig0            = ig0 + 1
408            ENDDO
409         ENDDO
410
411         pksurcp       = ppk(1,jjp1,l) / cpp
412         zplay(ig0,l)  = preff * pksurcp ** unskap
413!         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
414
415      ENDDO
416
417      if (flag_moyzon) then
418        zplaybar(:,:) = preff * (zpkbar(:,:)/cpp)**unskap
419      endif
420
421c   43.bis traceurs (tous intensifs) 
422c   ---------------
423
424      DO iq=1,nqtot
425         DO l=1,llm
426            zqfi(1,l,iq) = pq(1,1,l,iq)
427            ig0          = 2
428            DO j=2,jjm
429               DO i = 1, iim
430                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
431                  ig0             = ig0 + 1
432               ENDDO
433            ENDDO
434            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
435         ENDDO
436      ENDDO  ! boucle sur traceurs
437
438      if (flag_moyzon) then
439       DO iq=1,nqtot
440! RQ: REVOIR A QUOI CA SERT... ET VERIFIER...
441!       iiq=niadv(iq)
442! en fait, iiq=iq...
443! FIN RQ
444        tmpvar(:,:) = pq(:,1,:,iq)
445        call moyzon(llm,tmpvar,tmpvarbar)
446        zqfibar(1,:,iq) = tmpvarbar
447        ig0 = 2
448         do j = 2, jjm
449          tmpvar(:,:) = pq(:,j,:,iq)
450          call moyzon(llm,tmpvar,tmpvarbar)
451          zqfibar(ig0,:,iq) = tmpvarbar
452          ig0 = ig0+1
453          do i=2,iim
454            zqfibar(ig0,:,iq) = zqfibar(ig0-1,:,iq)
455            ig0 = ig0+1
456          enddo
457         enddo
458        tmpvar(:,:) = pq(:,jjp1,:,iq)
459        call moyzon(llm,tmpvar,tmpvarbar)
460        zqfibar(ngridmx,:,iq) = tmpvarbar
461       ENDDO ! of DO iq=1,nqtot
462      endif
463
464! DEBUG
465!     do ig0=1,ngridmx
466!       write(*,'(6(e13.5,1x))') zqfibar(ig0,1,10),zqfi(ig0,1,10),
467!    .                         zqfibar(ig0,llm/2,10),zqfi(ig0,llm/2,10),
468!    .                           zqfibar(ig0,llm,10),zqfi(ig0,llm,10)
469!     enddo
470!     stop
471
472!-----------------
473! RQ SL 13/10/10:
474! Ces calculs ne servent pas.
475! Si necessaire, decommenter ces variables et les calculs...
476!
477!   convergence dynamique pour les traceurs "EAU"
478! Earth-specific treatment of first 2 tracers (water)
479!      if (planet_type=="earth") then
480!       DO iq=1,2
481!        DO l=1,llm
482!           pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
483!           ig0          = 2
484!           DO j=2,jjm
485!              DO i = 1, iim
486!                 pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
487!                 ig0             = ig0 + 1
488!              ENDDO
489!           ENDDO
490!           pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
491!        ENDDO
492!       ENDDO
493!      endif ! of if (planet_type=="earth")
494!----------------
495
496c   Geopotentiel calcule par rapport a la surface locale:
497c   -----------------------------------------------------
498
499      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
500      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
501      DO l=1,llm
502         DO ig=1,ngridmx
503           zphi(ig,l)=zphi(ig,l)-zphis(ig)
504         ENDDO
505      ENDDO
506
507      if (flag_moyzon) then
508        tmpvar(:,1) = pphis(:,1)
509        call moyzon(1,tmpvar(:,1),tmpvarbar(1))
510        zphisbar(1) = tmpvarbar(1)
511        tmpvar(:,:) = pphi(:,1,:)
512        call moyzon(llm,tmpvar,tmpvarbar)
513        zphibar(1,:) = tmpvarbar
514        ig0 = 2
515         do j = 2, jjm
516          tmpvar(:,1) = pphis(:,j)
517          call moyzon(1,tmpvar(:,1),tmpvarbar(1))
518          zphisbar(ig0) = tmpvarbar(1)
519          tmpvar(:,:) = pphi(:,j,:)
520          call moyzon(llm,tmpvar,tmpvarbar)
521          zphibar(ig0,:) = tmpvarbar
522          ig0 = ig0+1
523          do i=2,iim
524            zphisbar(ig0)  = zphisbar(ig0-1)
525            zphibar(ig0,:) = zphibar(ig0-1,:)
526            ig0 = ig0+1
527          enddo
528         enddo
529        tmpvar(:,1) = pphis(:,jjp1)
530        call moyzon(1,tmpvar(:,1),tmpvarbar(1))
531        zphisbar(ngridmx) = tmpvarbar(1)
532        tmpvar(:,:) = pphi(:,jjp1,:)
533        call moyzon(llm,tmpvar,tmpvarbar)
534        zphibar(ngridmx,:) = tmpvarbar
535      endif
536
537c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
538c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux 
539c    de masse est calclue dans advtrac.F 
540c      DO l=1,llm
541c        pvervel(1,l)=pw(1,1,l) * g /apoln
542c        ig0=2
543c       DO j=2,jjm
544c           DO i = 1, iim
545c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
546c              ig0 = ig0 + 1
547c           ENDDO
548c       ENDDO
549c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
550c      ENDDO
551
552c
553c   45. champ u:
554c   ------------
555
556      DO 50 l=1,llm
557
558         DO 25 j=2,jjm
559            ig0 = 1+(j-2)*iim
560            zufi(ig0+1,l)= 0.5 * 
561     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
562!            pcvgu(ig0+1,l)= 0.5 *
563!     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
564            DO 10 i=2,iim
565               zufi(ig0+i,l)= 0.5 *
566     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
567!               pcvgu(ig0+i,l)= 0.5 *
568!     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
56910         CONTINUE
57025      CONTINUE
571
57250    CONTINUE
573
574
575c   46.champ v:
576c   -----------
577
578      DO l=1,llm
579         DO j=2,jjm
580            ig0=1+(j-2)*iim
581            DO i=1,iim
582               zvfi(ig0+i,l)= 0.5 *
583     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
584c               pcvgv(ig0+i,l)= 0.5 *
585c     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
586            ENDDO
587         ENDDO
588      ENDDO
589
590
591c   47. champs de vents aux pole nord   
592c   ------------------------------
593c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
594c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
595
596      DO l=1,llm
597
598         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
599         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
600         DO i=2,iim
601            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
602            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
603         ENDDO
604
605         DO i=1,iim
606            zcos(i)   = COS(rlonv(i))*z1(i)
607            zcosbis(i)= COS(rlonv(i))*z1bis(i)
608            zsin(i)   = SIN(rlonv(i))*z1(i)
609            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
610         ENDDO
611
612         zufi(1,l)  = SSUM(iim,zcos,1)/pi
613!         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
614         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
615!         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
616
617      ENDDO
618
619
620c   48. champs de vents aux pole sud:
621c   ---------------------------------
622c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
623c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
624
625      DO l=1,llm
626
627         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
628         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
629         DO i=2,iim
630            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
631            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
632         ENDDO
633
634         DO i=1,iim
635            zcos(i)    = COS(rlonv(i))*z1(i)
636            zcosbis(i) = COS(rlonv(i))*z1bis(i)
637            zsin(i)    = SIN(rlonv(i))*z1(i)
638            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
639         ENDDO
640
641         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
642!         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
643         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
644!         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
645
646      ENDDO
647c
648c On change de grille, dynamique vers physiq, pour le flux de masse verticale
649      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
650
651c-----------------------------------------------------------------------
652c   Appel de la physique:
653c   ---------------------
654
655! Appel de la physique: pose probleme quand on tourne
656! SANS physique, car physiq.F est dans le repertoire phy[]...
657! Il faut une cle CPP_PHYS
658
659! Le fait que les arguments de physiq soient differents selon les planetes
660! ne pose pas de probleme a priori.
661
662!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
663      zdt_split=dtphys/nsplit_phys
664      zdufic(:,:)=0.
665      zdvfic(:,:)=0.
666      zdtfic(:,:)=0.
667      zdqfic(:,:,:)=0.
668
669#ifdef CPP_PHYS
670
671      do isplit=1,nsplit_phys
672
673         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
674         debut_split=debut.and.isplit==1
675         lafin_split=lafin.and.isplit==nsplit_phys
676
677      if (planet_type.eq."earth") then
678         CALL physiq (ngridmx,
679     .             llm,
680     .             debut_split,
681     .             lafin_split,
682     .             jD_cur,
683     .             jH_cur_split,
684     .             zdt_split,
685     .             zplev,
686     .             zplay,
687     .             zphi,
688     .             zphis,
689     .             presnivs,
690     .             zufi,
691     .             zvfi,
692     .             ztfi,
693     .             zqfi,
694     .             flxwfi,
695     .             zdufi,
696     .             zdvfi,
697     .             zdtfi,
698     .             zdqfi,
699     .             zdpsrf,
700     .             pducov)
701
702      else if ( planet_type=="generic" ) then
703
704         CALL physiq (ngridmx,     !! ngrid
705     .             llm,            !! nlayer
706     .             nqtot,          !! nq
707     .             tname,          !! tracer names from dynamical core (given in infotrac)
708     .             debut_split,    !! firstcall
709     .             lafin_split,    !! lastcall
710     .             jD_cur,         !! pday. see leapfrog
711     .             jH_cur_split,   !! ptime "fraction of day"
712     .             zdt_split,      !! ptimestep
713     .             zplev,          !! pplev
714     .             zplay,          !! pplay
715     .             zphi,           !! pphi
716     .             zufi,           !! pu
717     .             zvfi,           !! pv
718     .             ztfi,           !! pt
719     .             zqfi,           !! pq
720     .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
721     .             zdufi,          !! pdu
722     .             zdvfi,          !! pdv
723     .             zdtfi,          !! pdt
724     .             zdqfi,          !! pdq
725     .             zdpsrf,         !! pdpsrf
726     .             tracerdyn)      !! tracerdyn <-- utilite ???
727
728      else if ( planet_type=="mars" ) then
729
730        CALL physiq (ngridmx,    ! ngrid
731     .             llm,          ! nlayer
732     .             nqtot,        ! nq
733     .             debut_split,  ! firstcall
734     .             lafin_split,  ! lastcall
735     .             jD_cur,       ! pday
736     .             jH_cur_split, ! ptime
737     .             zdt_split,    ! ptimestep
738     .             zplev,        ! pplev
739     .             zplay,        ! pplay
740     .             zphi,         ! pphi
741     .             zufi,         ! pu
742     .             zvfi,         ! pv
743     .             ztfi,         ! pt
744     .             zqfi,         ! pq
745     .             flxwfi,       ! pw
746     .             zdufi,        ! pdu
747     .             zdvfi,        ! pdv
748     .             zdtfi,        ! pdt
749     .             zdqfi,        ! pdq
750     .             zdpsrf,       ! pdpsrf
751     .             tracerdyn)    ! tracerdyn (somewhat obsolete)
752
753      else if ((planet_type=="titan").or.(planet_type=="venus")) then
754
755         CALL physiq (ngridmx,
756     .             llm,
757     .             nqtot,
758     .             debut_split,
759     .             lafin_split,
760     .             jD_cur,
761     .             jH_cur_split,
762     .             zdt_split,
763     .             zplev,
764     .             zplay,
765     .             zpk,
766     .             zphi,
767     .             zphis,
768     .             presnivs,
769     .             zufi,
770     .             zvfi,
771     .             ztfi,
772     .             zqfi,
773     .             flxwfi,
774     .             zdufi,
775     .             zdvfi,
776     .             zdtfi,
777     .             zdqfi,
778     .             zdpsrf)
779
780      else ! unknown "planet_type"
781
782        write(lunout,*) "calfis_p: error, unknown planet_type: ",
783     &                  trim(planet_type)
784        stop
785
786      endif ! planet_type
787
788         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
789         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
790         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
791         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
792
793         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
794         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
795         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
796         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
797
798      enddo ! of do isplit=1,nsplit_phys
799
800! ATTENTION...
801      if (flag_moyzon.and.(nsplit_phys.ne.1)) then
802       print*,"WARNING ! flag_moyzon + nsplit_phys"
803       print*,"zqfibar n'est pas implemente au cours des iterations"
804       print*,"Donc a revoir..."
805       stop
806      endif
807
808#endif
809! #endif of #ifdef CPP_PHYS
810
811      zdufi(:,:)=zdufic(:,:)/nsplit_phys
812      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
813      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
814      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
815
816
817500   CONTINUE
818
819c-----------------------------------------------------------------------
820c   transformation des tendances physiques en tendances dynamiques:
821c   ---------------------------------------------------------------
822
823c  tendance sur la pression :
824c  -----------------------------------
825
826      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
827c
828c   62. enthalpie potentielle
829c   ---------------------
830
831! ADAPTATION GCM POUR CP(T)
832      call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
833
834         DO i=1,iip1
835      pdhfi(i,1,:) = (zteta(1,:) - pteta(i,1,:))/dtphys
836      pdhfi(i,jjp1,:) = (zteta(ngridmx,:) - pteta(i,jjp1,:))/dtphys
837         ENDDO
838
839         DO j=2,jjm
840            ig0=1+(j-2)*iim
841            DO i=1,iim
842      pdhfi(i,j,:) = (zteta(ig0+i,:) - pteta(i,j,:))/dtphys
843            ENDDO
844               pdhfi(iip1,j,:) =  pdhfi(1,j,:)
845         ENDDO
846
847
848c   62. humidite specifique
849c   ---------------------
850! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
851!      DO iq=1,nqtot
852!         DO l=1,llm
853!            DO i=1,iip1
854!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
855!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
856!            ENDDO
857!            DO j=2,jjm
858!               ig0=1+(j-2)*iim
859!               DO i=1,iim
860!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
861!               ENDDO
862!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
863!            ENDDO
864!         ENDDO
865!      ENDDO
866
867c   63. traceurs (tous en intensifs)
868c   ------------
869C     initialisation des tendances
870      pdqfi(:,:,:,:)=0.
871C
872       DO iq=1,nqtot
873         iiq=niadv(iq)
874         DO l=1,llm
875            DO i=1,iip1
876               pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
877               pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
878            ENDDO
879            DO j=2,jjm
880               ig0=1+(j-2)*iim
881               DO i=1,iim
882                  pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
883               ENDDO
884               pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
885            ENDDO
886         ENDDO
887       ENDDO
888
889c   65. champ u:
890c   ------------
891
892      DO l=1,llm
893
894         DO i=1,iip1
895            pdufi(i,1,l)    = 0.
896            pdufi(i,jjp1,l) = 0.
897         ENDDO
898
899         DO j=2,jjm
900            ig0=1+(j-2)*iim
901            DO i=1,iim-1
902               pdufi(i,j,l)=
903     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
904            ENDDO
905            pdufi(iim,j,l)=
906     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
907            pdufi(iip1,j,l)=pdufi(1,j,l)
908         ENDDO
909
910      ENDDO
911
912
913c   67. champ v:
914c   ------------
915
916      DO l=1,llm
917
918         DO j=2,jjm-1
919            ig0=1+(j-2)*iim
920            DO i=1,iim
921               pdvfi(i,j,l)=
922     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
923            ENDDO
924            pdvfi(iip1,j,l) = pdvfi(1,j,l)
925         ENDDO
926      ENDDO
927
928
929c   68. champ v pres des poles:
930c   ---------------------------
931c      v = U * cos(long) + V * SIN(long)
932
933      DO l=1,llm
934
935         DO i=1,iim
936            pdvfi(i,1,l)=
937     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
938            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
939     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
940            pdvfi(i,1,l)=
941     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
942            pdvfi(i,jjm,l)=
943     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
944          ENDDO
945
946         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
947         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
948
949      ENDDO
950
951c-----------------------------------------------------------------------
952
953700   CONTINUE
954 
955      firstcal = .FALSE.
956
957      RETURN
958      END
Note: See TracBrowser for help on using the repository browser.