source: trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/calfis.F @ 1477

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