source: trunk/LMDZ.MARS/libf/dynphy_lonlat/calfis.F @ 2414

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