source: trunk/LMDZ.TITAN/libf/dynphy_lonlat/calfis.F @ 3094

Last change on this file since 3094 was 1576, checked in by emillour, 8 years ago

All models:
Clean up the dynamics/physics interface to converge with LMDZ5;
get rid of tracerdyn parameter (which was supposed to send information
about wether tracers should be advected or not in the Mars and Generic
models).
EM

File size: 15.5 KB
Line 
1      SUBROUTINE calfis(nq, lafin, rdayvrai,rday_ecri, heure,
2     $            pucov,pvcov,pteta,pq,pmasse,pps,pp,ppk,pphis,pphi,
3     $            pducov,pdvcov,pdteta,pdq,pw,
4     $            pdufi,pdvfi,pdhfi,pdqfi,pdpsfi )
5c
6c    Auteur :  P. Le Van, F. Hourdin
7c   .........
8      USE infotrac, ONLY: tname, nqtot
9      USE comvert_mod, ONLY: preff
10      USE comconst_mod, ONLY: dtphys,kappa,cpp,pi
11      USE physiq_mod, ONLY: physiq
12
13      IMPLICIT NONE
14c=======================================================================
15c
16c   1. rearrangement des tableaux et transformation
17c      variables dynamiques  >  variables physiques
18c   2. calcul des termes physiques
19c   3. retransformation des tendances physiques en tendances dynamiques
20c
21c   remarques:
22c   ----------
23c
24c    - les vents sont donnes dans la physique par leurs composantes
25c      naturelles.
26c    - la variable thermodynamique de la physique est une variable
27c      intensive :   T
28c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
29c    - les deux seules variables dependant de la geometrie necessaires
30c      pour la physique sont la latitude pour le rayonnement et
31c      l'aire de la maille quand on veut integrer une grandeur
32c      horizontalement.
33c    - les points de la physique sont les points scalaires de la
34c      la dynamique; numerotation:
35c          1 pour le pole nord
36c          (jjm-1)*iim pour l'interieur du domaine
37c          ngridmx pour le pole sud
38c      ---> ngridmx=2+(jjm-1)*iim
39c
40c     Input :
41c     -------
42c       ecritphy        frequence d'ecriture (en jours)de histphy
43c       pucov           covariant zonal velocity
44c       pvcov           covariant meridional velocity
45c       pteta           potential temperature
46c       pps             surface pressure
47c       pmasse          masse d'air dans chaque maille
48c       pts             surface temperature  (K)
49c       pw              flux vertical (kg/s)
50c
51c    Output :
52c    --------
53c        pdufi          tendency for the natural zonal velocity (ms-1)
54c        pdvfi          tendency for the natural meridional velocity
55c        pdhfi          tendency for the potential temperature
56c        pdtsfi         tendency for the surface temperature
57c
58c=======================================================================
59c
60c-----------------------------------------------------------------------
61c
62c    0.  Declarations :
63c    ------------------
64
65#include "dimensions.h"
66#include "paramet.h"
67
68      INTEGER ngridmx,nq
69      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
70
71#include "comgeom2.h"
72!#include "control.h"
73
74!#include "advtrac.h"
75!! this is to get tnom (tracers name)
76
77c    Arguments :
78c    -----------
79      LOGICAL  lafin
80      REAL heure
81
82      REAL pvcov(iip1,jjm,llm)
83      REAL pucov(iip1,jjp1,llm)
84      REAL pteta(iip1,jjp1,llm)
85      REAL pmasse(iip1,jjp1,llm)
86      REAL pq(iip1,jjp1,llm,nqtot)
87      REAL pphis(iip1,jjp1)
88      REAL pphi(iip1,jjp1,llm)
89c
90      REAL pdvcov(iip1,jjm,llm)
91      REAL pducov(iip1,jjp1,llm)
92      REAL pdteta(iip1,jjp1,llm)
93      REAL pdq(iip1,jjp1,llm,nqtot)
94c
95      REAL pw(iip1,jjp1,llm)
96c
97      REAL pps(iip1,jjp1)
98      REAL pp(iip1,jjp1,llmp1)
99      REAL ppk(iip1,jjp1,llm)
100c
101      REAL pdvfi(iip1,jjm,llm)
102      REAL pdufi(iip1,jjp1,llm)
103      REAL pdhfi(iip1,jjp1,llm)
104      REAL pdqfi(iip1,jjp1,llm,nqtot)
105      REAL pdpsfi(iip1,jjp1)
106
107c    Local variables :
108c    -----------------
109
110      INTEGER i,j,l,ig0,ig,iq
111      REAL zpsrf(ngridmx)
112      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
113      REAL zphi(ngridmx,llm),zphis(ngridmx)
114c
115      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
116      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
117c
118!      REAL zvervel(ngridmx,llm)
119      REAL flxwfi(ngridmx,llm) ! vertical mass flux (kg/s) on physics grid
120c
121      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
122      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
123      REAL zdpsrf(ngridmx)
124c
125      REAL zsin(iim),zcos(iim),z1(iim)
126      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
127      REAL unskap, pksurcp
128c
129     
130      EXTERNAL gr_dyn_fi,gr_fi_dyn
131      REAL SSUM
132      EXTERNAL SSUM
133
134      REAL latfi(ngridmx),lonfi(ngridmx)
135      REAL airefi(ngridmx)
136      SAVE latfi, lonfi, airefi
137
138      LOGICAL firstcal, debut
139      DATA firstcal/.true./
140      SAVE firstcal,debut
141      REAL rdayvrai,rday_ecri
142c
143c-----------------------------------------------------------------------
144c
145c    1. Initialisations :
146c    --------------------
147c
148
149
150      IF (ngridmx.NE.2+(jjm-1)*iim) THEN
151         PRINT*,'STOP dans calfis'
152         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
153         PRINT*,'  ngridmx  jjm   iim   '
154         PRINT*,ngridmx,jjm,iim
155         STOP
156      ENDIF
157
158c-----------------------------------------------------------------------
159c   latitude, longitude et aires des mailles pour la physique:
160c   ----------------------------------------------------------
161
162c
163      IF ( firstcal )  THEN
164          debut = .TRUE.
165      ELSE
166          debut = .FALSE.
167      ENDIF
168
169c
170!      IF (firstcal) THEN
171!         latfi(1)=rlatu(1)
172!         lonfi(1)=0.
173!         DO j=2,jjm
174!            DO i=1,iim
175!               latfi((j-2)*iim+1+i)= rlatu(j)
176!               lonfi((j-2)*iim+1+i)= rlonv(i)
177!            ENDDO
178!         ENDDO
179!         latfi(ngridmx)= rlatu(jjp1)
180!         lonfi(ngridmx)= 0.
181!
182!         ! build airefi(), mesh area on physics grid
183!         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
184!         ! Poles are single points on physics grid
185!         airefi(1)=airefi(1)*iim
186!         airefi(ngridmx)=airefi(ngridmx)*iim
187!
188!         CALL inifis(ngridmx,llm,day_ini,daysec,dtphys,
189!     .                latfi,lonfi,airefi,rad,g,r,cpp)
190!      ENDIF
191
192c
193c-----------------------------------------------------------------------
194c   40. transformation des variables dynamiques en variables physiques:
195c   ---------------------------------------------------------------
196
197c   41. pressions au sol (en Pascals)
198c   ----------------------------------
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
212c   42. pression intercouches :
213c
214c   -----------------------------------------------------------------
215c     .... zplev  definis aux (llm +1) interfaces des couches  ....
216c     .... zplay  definis aux (  llm )    milieux des couches  ....
217c   -----------------------------------------------------------------
218
219c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
220c
221       unskap   = 1./ kappa
222c
223      DO l = 1, llmp1
224        zplev( 1,l ) = pp(1,1,l)
225        ig0 = 2
226          DO j = 2, jjm
227             DO i =1, iim
228              zplev( ig0,l ) = pp(i,j,l)
229              ig0 = ig0 +1
230             ENDDO
231          ENDDO
232        zplev( ngridmx,l ) = pp(1,jjp1,l)
233      ENDDO
234c
235c
236
237c   43. temperature naturelle (en K) et pressions milieux couches .
238c   ---------------------------------------------------------------
239
240      DO l=1,llm
241
242         pksurcp     =  ppk(1,1,l) / cpp
243         zplay(1,l)  =  preff * pksurcp ** unskap
244         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
245         ig0         =  2
246
247         DO j = 2, jjm
248            DO i = 1, iim
249              pksurcp        = ppk(i,j,l) / cpp
250              zplay(ig0,l)   = preff * pksurcp ** unskap
251              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
252              ig0            = ig0 + 1
253            ENDDO
254         ENDDO
255
256         pksurcp       = ppk(1,jjp1,l) / cpp
257         zplay(ig0,l)  = preff * pksurcp ** unskap
258         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
259
260      ENDDO
261
262      DO l=1, llm
263        DO ig=1,ngridmx
264             if (ztfi(ig,l).lt.15) then
265                  write(*,*) 'New Temperature below 15 K !!! '
266                      write(*,*) 'Stop in calfis.F '
267                      write(*,*) 'ig=', ig, ' l=', l
268                      write(*,*) 'ztfi(ig,l)=',ztfi(ig,l)
269                  stop
270             end if
271        ENDDO
272      ENDDO
273
274
275
276c   43.bis Taceurs (en kg/kg)
277c   --------------------------
278      DO iq=1,nqtot
279         DO l=1,llm
280            zqfi(1,l,iq) = pq(1,1,l,iq)
281            ig0          = 2
282            DO j=2,jjm
283               DO i = 1, iim
284                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
285                  ig0             = ig0 + 1
286               ENDDO
287            ENDDO
288            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
289         ENDDO
290      ENDDO
291
292c   Geopotentiel calcule par rapport a la surface locale:
293c   -----------------------------------------------------
294
295      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
296      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
297      DO l=1,llm
298         DO ig=1,ngridmx
299            zphi(ig,l)=zphi(ig,l)-zphis(ig)
300         ENDDO
301      ENDDO
302
303c   Calcul de la vitesse  verticale (m/s) pour diagnostique
304c   -------------------------------
305c     pw est en kg/s
306c On interpole "lineairement" la temperature entre les couches(FF,10/95)
307
308!      DO ig=1,ngridmx
309!         zvervel(ig,1)=0.
310!      END DO
311!      DO l=2,llm
312!        zvervel(1,l)=(pw(1,1,l)/apoln)
313!     &  * r *0.5*(ztfi(1,l)+ztfi(1,l-1)) /zplev(1,l)             
314!        ig0=2
315!       DO j=2,jjm
316!           DO i = 1, iim
317!              zvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
318!     &        * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
319!              ig0 = ig0 + 1
320!           ENDDO
321!       ENDDO
322!        zvervel(ig0,l)=(pw(1,jjp1,l)/apols)
323!     &  * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l)             
324!      ENDDO
325
326c    .........  Reindexation : calcul de zvervel au MILIEU des couches
327!       DO l=1,llm-1
328!             DO ig=1,ngridmx
329!                    zvervel(ig,l) = 0.5*(zvervel(ig,l)+zvervel(ig,l+1))
330!          END DO
331!       END DO
332c      (dans la couche llm, on garde la valeur à la limite inférieure llm)
333
334! vertical mass flux
335      ! tranfer values from dynamics grid to physics grid:
336      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pw,flxwfi)
337      ! but mass flux is an extensive variable, so take the sum at the poles
338      DO l=1,llm
339        flxwfi(1,l)=sum(pw(1:iim,1,l))
340        flxwfi(ngridmx,l)=sum(pw(1:iim,jjp1,l))
341      ENDDO
342
343c   45. champ u:
344c   ------------
345
346      DO l=1,llm
347         DO j=2,jjm
348            ig0 = 1+(j-2)*iim
349            zufi(ig0+1,l)= 0.5 *
350     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
351            DO i=2,iim
352               zufi(ig0+i,l)= 0.5 *
353     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
354            ENDDO
355        ENDDO
356      ENDDO
357
358
359c   46.champ v:
360c   -----------
361
362      DO l=1,llm
363         DO j=2,jjm
364            ig0=1+(j-2)*iim
365            DO i=1,iim
366               zvfi(ig0+i,l)= 0.5 *
367     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
368            ENDDO
369         ENDDO
370      ENDDO
371
372
373c   47. champs de vents aux pole nord   
374c   ------------------------------
375c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
376c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
377
378      DO l=1,llm
379
380         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
381         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
382         DO i=2,iim
383            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
384            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
385         ENDDO
386
387         DO i=1,iim
388            zcos(i)   = COS(rlonv(i))*z1(i)
389            zcosbis(i)= COS(rlonv(i))*z1bis(i)
390            zsin(i)   = SIN(rlonv(i))*z1(i)
391            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
392         ENDDO
393
394         zufi(1,l)  = SSUM(iim,zcos,1)/pi
395         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
396
397      ENDDO
398
399
400c   48. champs de vents aux pole sud:
401c   ---------------------------------
402c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
403c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
404
405      DO l=1,llm
406
407         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
408         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
409         DO i=2,iim
410            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
411            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
412      ENDDO
413
414         DO i=1,iim
415            zcos(i)    = COS(rlonv(i))*z1(i)
416            zcosbis(i) = COS(rlonv(i))*z1bis(i)
417            zsin(i)    = SIN(rlonv(i))*z1(i)
418            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
419      ENDDO
420
421         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
422         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
423
424      ENDDO
425
426c-----------------------------------------------------------------------
427c   Appel de la physique:
428c   ---------------------
429
430
431      CALL physiq (ngridmx,llm,nq,
432     .     tname,
433     ,     debut,lafin,
434     ,     rday_ecri,heure,dtphys,
435     ,     zplev,zplay,zphi,
436     ,     zufi, zvfi,ztfi, zqfi, 
437!     ,     zvervel,
438     ,     flxwfi,
439C - sorties
440     s     zdufi, zdvfi, zdtfi, zdqfi,zdpsrf)
441
442
443c-----------------------------------------------------------------------
444c   transformation des tendances physiques en tendances dynamiques:
445c   ---------------------------------------------------------------
446
447c  tendance sur la pression :
448c  -----------------------------------
449
450      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
451c
452ccc     CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
453
454c   62. enthalpie potentielle
455c   ---------------------
456
457      DO l=1,llm
458
459         DO i=1,iip1
460          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
461          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
462         ENDDO
463
464         DO j=2,jjm
465            ig0=1+(j-2)*iim
466            DO i=1,iim
467               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
468            ENDDO
469               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
470         ENDDO
471
472      ENDDO
473
474
475c   62. traceurs
476c   ---------------------
477
478      DO iq=1,nqtot
479         DO l=1,llm
480            DO i=1,iip1
481               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
482               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
483            ENDDO
484            DO j=2,jjm
485               ig0=1+(j-2)*iim
486               DO i=1,iim
487                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
488               ENDDO
489               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
490            ENDDO
491         ENDDO
492      ENDDO
493
494c   65. champ u:
495c   ------------
496
497      DO l=1,llm
498
499         DO i=1,iip1
500            pdufi(i,1,l)    = 0.
501            pdufi(i,jjp1,l) = 0.
502         ENDDO
503
504         DO j=2,jjm
505            ig0=1+(j-2)*iim
506            DO i=1,iim-1
507               pdufi(i,j,l)=
508     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
509            ENDDO
510            pdufi(iim,j,l)=
511     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
512            pdufi(iip1,j,l)=pdufi(1,j,l)
513         ENDDO
514
515      ENDDO
516
517
518c   67. champ v:
519c   ------------
520
521      DO l=1,llm
522
523         DO j=2,jjm-1
524            ig0=1+(j-2)*iim
525            DO i=1,iim
526               pdvfi(i,j,l)=
527     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
528            ENDDO
529            pdvfi(iip1,j,l) = pdvfi(1,j,l)
530         ENDDO
531      ENDDO
532
533
534c   68. champ v pres des poles:
535c   ---------------------------
536c      v = U * cos(long) + V * SIN(long)
537
538      DO l=1,llm
539
540         DO i=1,iim
541            pdvfi(i,1,l)=
542     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
543            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
544     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
545            pdvfi(i,1,l)=
546     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
547            pdvfi(i,jjm,l)=
548     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
549          ENDDO
550
551         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
552         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
553
554      ENDDO
555
556c-----------------------------------------------------------------------
557
558700   CONTINUE
559
560      firstcal = .FALSE.
561
562      RETURN
563      END
Note: See TracBrowser for help on using the repository browser.