source: trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/calfis.F @ 1772

Last change on this file since 1772 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

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