source: trunk/libf/dyn3d/leapfrog.F @ 110

Last change on this file since 110 was 110, checked in by slebonnois, 14 years ago

SL: corrections de bugs suite a compilations venus et titan de la version 109.

File size: 27.3 KB
Line 
1!
2! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
3!
4c
5c
6      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,
7     &                    time_0)
8
9
10cIM : pour sortir les param. du modele dans un fis. netcdf 110106
11#ifdef CPP_IOIPSL
12      use IOIPSL
13#endif
14      USE infotrac
15      USE guide_mod, ONLY : guide_main
16      USE write_field
17      USE control_mod
18      IMPLICIT NONE
19
20c      ......   Version  du 10/01/98    ..........
21
22c             avec  coordonnees  verticales hybrides
23c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
24
25c=======================================================================
26c
27c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
28c   -------
29c
30c   Objet:
31c   ------
32c
33c   GCM LMD nouvelle grille
34c
35c=======================================================================
36c
37c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
38c      et possibilite d'appeler une fonction f(y)  a derivee tangente
39c      hyperbolique a la  place de la fonction a derivee sinusoidale.
40
41c  ... Possibilite de choisir le shema pour l'advection de
42c        q  , en modifiant iadv dans traceur.def  (10/02) .
43c
44c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
45c      Pour Van-Leer iadv=10
46c
47c-----------------------------------------------------------------------
48c   Declarations:
49c   -------------
50
51#include "dimensions.h"
52#include "paramet.h"
53#include "comconst.h"
54#include "comdissnew.h"
55#include "comvert.h"
56#include "comgeom.h"
57#include "logic.h"
58#include "temps.h"
59#include "ener.h"
60#include "description.h"
61#include "serre.h"
62!#include "com_io_dyn.h"
63#include "iniprint.h"
64#include "academic.h"
65
66! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
67! #include "clesphys.h"
68
69      real zqmin,zqmax
70      INTEGER nbetatmoy, nbetatdem,nbetat
71
72c   variables dynamiques
73      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
74      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
75      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
76      REAL ps(ip1jmp1)                       ! pression  au sol
77      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
78      REAL pks(ip1jmp1)                      ! exner au  sol
79      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
80      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
81      REAL masse(ip1jmp1,llm)                ! masse d'air
82      REAL phis(ip1jmp1)                     ! geopotentiel au sol
83      REAL phi(ip1jmp1,llm)                  ! geopotentiel
84      REAL w(ip1jmp1,llm)                    ! vitesse verticale
85! ADAPTATION GCM POUR CP(T)
86      REAL temp(ip1jmp1,llm)                 ! temperature 
87      REAL tsurpk(ip1jmp1,llm)               ! cpp*T/pk 
88
89c variables dynamiques intermediaire pour le transport
90      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
91
92c   variables dynamiques au pas -1
93      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
94      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
95      REAL massem1(ip1jmp1,llm)
96
97c   tendances dynamiques en */s
98      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
99      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
100
101c   tendances de la dissipation en */s
102      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
103      REAL dtetadis(ip1jmp1,llm)
104
105c   tendances de la couche superieure */s
106      REAL dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
107      REAL dtetatop(ip1jmp1,llm)
108      REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
109
110c   tendances physiques */s
111      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
112      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
113
114c   variables pour le fichier histoire
115      REAL dtav      ! intervalle de temps elementaire
116
117      REAL tppn(iim),tpps(iim),tpn,tps
118c
119      INTEGER itau,itaufinp1,iav
120!      INTEGER  iday ! jour julien
121      REAL       time
122
123      REAL  SSUM
124      REAL time_0 , finvmaold(ip1jmp1,llm)
125
126cym      LOGICAL  lafin
127      LOGICAL :: lafin=.false.
128      INTEGER ij,iq,l
129      INTEGER ik
130
131      real time_step, t_wrt, t_ops
132
133!      REAL rdayvrai,rdaym_ini
134! jD_cur: jour julien courant
135! jH_cur: heure julienne courante
136      REAL :: jD_cur, jH_cur
137      INTEGER :: an, mois, jour
138      REAL :: secondes
139
140      LOGICAL first,callinigrads
141cIM : pour sortir les param. du modele dans un fis. netcdf 110106
142      save first
143      data first/.true./
144      real dt_cum
145      character*10 infile
146      integer zan, tau0, thoriid
147      integer nid_ctesGCM
148      save nid_ctesGCM
149      real degres
150      real rlong(iip1), rlatg(jjp1)
151      real zx_tmp_2d(iip1,jjp1)
152      integer ndex2d(iip1*jjp1)
153      logical ok_sync
154      parameter (ok_sync = .true.)
155
156      data callinigrads/.true./
157      character*10 string10
158
159      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
160      REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
161
162c+jld variables test conservation energie
163      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
164C     Tendance de la temp. potentiel d (theta)/ d t due a la
165C     tansformation d'energie cinetique en energie thermique
166C     cree par la dissipation
167      REAL dtetaecdt(ip1jmp1,llm)
168      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
169      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
170      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
171      CHARACTER*15 ztit
172!IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
173!IM   SAVE      ip_ebil_dyn
174!IM   DATA      ip_ebil_dyn/0/
175c-jld
176
177      integer :: itau_w ! for write_paramLMDZ_dyn.h
178
179      character*80 dynhist_file, dynhistave_file
180      character(len=20) :: modname
181      character*80 abort_message
182
183      logical dissip_conservative
184      save dissip_conservative
185      data dissip_conservative/.true./
186
187      INTEGER testita
188      PARAMETER (testita = 9)
189
190      logical , parameter :: flag_verif = .false.
191     
192      ! for CP(T)
193      real :: dtec
194      real,external :: cpdet
195      real :: ztetaec(ip1jmp1,llm)
196
197c dummy: sinon cette routine n'est jamais compilee...
198      if(1.eq.0) then
199        CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
200      endif
201
202      itaufin   = nday*day_step
203      if (less1day) then
204c MODIF VENUS: to run less than one day:
205        itaufin   = int(fractday*day_step)
206      endif
207      itaufinp1 = itaufin +1
208      modname="leapfrog"
209     
210c INITIALISATIONS
211        dudis(:,:)   =0.
212        dvdis(:,:)   =0.
213        dtetadis(:,:)=0.
214        dutop(:,:)   =0.
215        dvtop(:,:)   =0.
216        dtetatop(:,:)=0.
217        dqtop(:,:,:) =0.
218        dptop(:)     =0.
219        dufi(:,:)   =0.
220        dvfi(:,:)   =0.
221        dtetafi(:,:)=0.
222        dqfi(:,:,:) =0.
223        dpfi(:)     =0.
224
225      itau = 0
226c      iday = day_ini+itau/day_step
227c      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
228c         IF(time.GT.1.) THEN
229c          time = time-1.
230c          iday = iday+1
231c         ENDIF
232
233
234c-----------------------------------------------------------------------
235c   On initialise la pression et la fonction d'Exner :
236c   --------------------------------------------------
237
238      dq(:,:,:)=0.
239      CALL pression ( ip1jmp1, ap, bp, ps, p       )
240      if (planet_type.eq."earth") then
241        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
242      else
243        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
244      endif
245
246c------------------
247c TEST PK MONOTONE
248c------------------
249      write(*,*) "Test PK"
250      do ij=1,ip1jmp1
251        do l=2,llm
252          if(pk(ij,l).gt.pk(ij,l-1)) then
253c           write(*,*) ij,l,pk(ij,l)
254            abort_message = 'PK non strictement decroissante'
255            call abort_gcm(modname,abort_message,1)
256c           write(*,*) "ATTENTION, Test PK deconnecté..."
257          endif
258        enddo
259      enddo
260      write(*,*) "Fin Test PK"
261c     stop
262c------------------
263
264c-----------------------------------------------------------------------
265c   Debut de l'integration temporelle:
266c   ----------------------------------
267
268   1  CONTINUE
269
270      jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
271      jH_cur = jH_ref +                                                 &
272     &          (itau * dtvr / daysec - int(itau * dtvr / daysec))
273
274
275#ifdef CPP_IOIPSL
276      if (ok_guide) then
277        call guide_main(itau,ucov,vcov,teta,q,masse,ps)
278      endif
279#endif
280
281
282c
283c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
284c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
285c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
286c     ENDIF
287c
288
289! Save fields obtained at previous time step as '...m1'
290      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
291      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
292      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
293      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
294      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
295
296      forward = .TRUE.
297      leapf   = .FALSE.
298      dt      =  dtvr
299
300c   ...    P.Le Van .26/04/94  ....
301
302      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
303      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
304
305   2  CONTINUE
306
307c-----------------------------------------------------------------------
308
309c   date:
310c   -----
311
312
313c   gestion des appels de la physique et des dissipations:
314c   ------------------------------------------------------
315c
316c   ...    P.Le Van  ( 6/02/95 )  ....
317
318      apphys = .FALSE.
319      statcl = .FALSE.
320      conser = .FALSE.
321      apdiss = .FALSE.
322
323      IF( purmats ) THEN
324      ! Purely Matsuno time stepping
325         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
326         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
327         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
328     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
329      ELSE
330      ! Leapfrog/Matsuno time stepping
331         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
332         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
333         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
334      END IF
335
336! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
337!          supress dissipation step
338      if (llm.eq.1) then
339        apdiss=.false.
340      endif
341
342c-----------------------------------------------------------------------
343c   calcul des tendances dynamiques:
344c   --------------------------------
345
346! ADAPTATION GCM POUR CP(T)
347      call tpot2t(ijp1llm,teta,temp,pk)
348      tsurpk = cpp*temp/pk
349      CALL geopot  ( ip1jmp1, tsurpk  , pk , pks,  phis  , phi   )
350
351      time = jD_cur + jH_cur
352      CALL caldyn
353     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,tsurpk,phis ,
354     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
355
356
357c-----------------------------------------------------------------------
358c   calcul des tendances advection des traceurs (dont l'humidite)
359c   -------------------------------------------------------------
360
361      IF( forward. OR . leapf )  THEN
362
363         CALL caladvtrac(q,pbaru,pbarv,
364     *        p, masse, dq,  teta,
365     .        flxw, pk)
366         
367         IF (offline) THEN
368Cmaf stokage du flux de masse pour traceurs OFF-LINE
369
370#ifdef CPP_IOIPSL
371           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
372     .   dtvr, itau)
373#endif
374
375
376         ENDIF ! of IF (offline)
377c
378      ENDIF ! of IF( forward. OR . leapf )
379
380
381c-----------------------------------------------------------------------
382c   integrations dynamique et traceurs:
383c   ----------------------------------
384
385
386       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
387     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
388     $              finvmaold                                    )
389
390
391c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
392c
393c-----------------------------------------------------------------------
394c   calcul des tendances physiques:
395c   -------------------------------
396c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
397c
398       IF( purmats )  THEN
399          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
400       ELSE
401          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
402       ENDIF
403c
404c
405       IF( apphys )  THEN
406c
407c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
408c
409
410         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
411         if (planet_type.eq."earth") then
412           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
413         else
414           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
415         endif
416
417!           rdaym_ini  = itau * dtvr / daysec
418!           rdayvrai   = rdaym_ini  + day_ini
419           jD_cur = jD_ref + day_ini - day_ref
420     $        + int (itau * dtvr / daysec)
421           jH_cur = jH_ref +                                            &
422     &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
423!         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
424!         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
425!         write(lunout,*)'current date = ',an, mois, jour, secondes
426
427c rajout debug
428c       lafin = .true.
429
430
431c   Interface avec les routines de phylmd (phymars ... )
432c   -----------------------------------------------------
433
434c+jld
435
436c  Diagnostique de conservation de l'énergie : initialisation
437         IF (ip_ebil_dyn.ge.1 ) THEN
438          ztit='bil dyn'
439! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
440           IF (planet_type.eq."earth") THEN
441            CALL diagedyn(ztit,2,1,1,dtphys
442     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
443           ENDIF
444         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
445c-jld
446#ifdef CPP_IOIPSL
447cIM : pour sortir les param. du modele dans un fis. netcdf 110106
448         IF (first) THEN
449#include "ini_paramLMDZ_dyn.h"
450         ENDIF
451c
452#include "write_paramLMDZ_dyn.h"
453c
454#endif
455! #endif of #ifdef CPP_IOIPSL
456
457         CALL calfis( lafin , jD_cur, jH_cur,
458     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
459     $               du,dv,dteta,dq,
460     $               flxw,
461     $               dufi,dvfi,dtetafi,dqfi,dpfi  )
462
463c      ajout des tendances physiques:
464c      ------------------------------
465          CALL addfi( dtphys, leapf, forward   ,
466     $                  ucov, vcov, teta , q   ,ps ,
467     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
468
469c      Couche superieure :
470c      -------------------
471         IF (ok_strato) THEN
472           CALL top_bound( vcov,ucov,teta,phi,masse,
473     $                     dutop,dvtop,dtetatop)
474c dqtop=0, dptop=0
475           CALL addfi( dtphys, leapf, forward   ,
476     $                  ucov, vcov, teta , q   ,ps ,
477     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
478         ENDIF
479
480c  Diagnostique de conservation de l'énergie : difference
481         IF (ip_ebil_dyn.ge.1 ) THEN
482          ztit='bil phys'
483          IF (planet_type.eq."earth") THEN
484           CALL diagedyn(ztit,2,1,1,dtphys
485     &     , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
486          ENDIF
487         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
488
489       ENDIF ! of IF( apphys )
490
491      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
492!   Academic case : Simple friction and Newtonan relaxation
493!   -------------------------------------------------------
494        DO l=1,llm   
495          DO ij=1,ip1jmp1
496           teta(ij,l)=teta(ij,l)-dtvr*
497     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
498          ENDDO
499        ENDDO ! of DO l=1,llm
500          call friction(ucov,vcov,dtvr)
501   
502       if (planet_type.eq."giant") then
503          ! Intrinsic heat flux
504          ! Aymeric -- for giant planets
505          if (ihf .gt. 1.e-6) then
506          !print *, '**** INTRINSIC HEAT FLUX ****', ihf
507            DO ij=1,ip1jmp1
508              teta(ij,1) = teta(ij,1)
509     &        + dtvr * aire(ij) * ihf / cpp / masse(ij,1)
510            ENDDO
511          !print *, '**** d teta '
512          !print *, dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
513          endif
514       endif
515
516   
517        ! Sponge layer (if any)
518        IF (ok_strato) THEN
519          CALL top_bound(vcov,ucov,teta,phi,
520     $                   masse,dutop,dvtop,dtetatop)
521c dqtop=0, dptop=0
522          CALL addfi( dtvr, leapf, forward   ,
523     $                  ucov, vcov, teta , q   ,ps ,
524     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
525        ENDIF ! of IF (ok_strato)
526      ENDIF ! of IF (iflag_phys.EQ.2)
527
528
529c-jld
530
531        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
532        if (planet_type.eq."earth") then
533          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
534        else
535          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
536        endif
537
538
539c-----------------------------------------------------------------------
540c   dissipation horizontale et verticale  des petites echelles:
541c   ----------------------------------------------------------
542
543      IF(apdiss) THEN
544
545
546c   calcul de l'energie cinetique avant dissipation
547        call covcont(llm,ucov,vcov,ucont,vcont)
548        call enercin(vcov,ucov,vcont,ucont,ecin0)
549
550c   dissipation
551! ADAPTATION GCM POUR CP(T)
552        call tpot2t(ijp1llm,teta,temp,pk)
553
554        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
555        ucov=ucov+dudis
556        vcov=vcov+dvdis
557        dudis=dudis/dtdiss   ! passage en (m/s)/s
558        dvdis=dvdis/dtdiss   ! passage en (m/s)/s
559
560c------------------------------------------------------------------------
561        if (dissip_conservative) then
562C       On rajoute la tendance due a la transform. Ec -> E therm. cree
563C       lors de la dissipation
564            call covcont(llm,ucov,vcov,ucont,vcont)
565            call enercin(vcov,ucov,vcont,ucont,ecin)
566! ADAPTATION GCM POUR CP(T)
567            do ij=1,ip1jmp1
568              do l=1,llm
569                dtec = (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
570                temp(ij,l) = temp(ij,l) + dtec
571              enddo
572            enddo
573            call t2tpot(ijp1llm,temp,ztetaec,pk)
574            dtetaecdt=ztetaec-teta
575            dtetadis=dtetadis+dtetaecdt
576        endif
577        teta=teta+dtetadis
578        dtetadis=dtetadis/dtdiss   ! passage en K/s
579c------------------------------------------------------------------------
580
581
582c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
583c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
584c
585
586        DO l  =  1, llm
587          DO ij =  1,iim
588           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
589           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
590          ENDDO
591           tpn  = SSUM(iim,tppn,1)/apoln
592           tps  = SSUM(iim,tpps,1)/apols
593
594          DO ij = 1, iip1
595           teta(  ij    ,l) = tpn
596           teta(ij+ip1jm,l) = tps
597          ENDDO
598        ENDDO
599
600        DO ij =  1,iim
601          tppn(ij)  = aire(  ij    ) * ps (  ij    )
602          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
603        ENDDO
604          tpn  = SSUM(iim,tppn,1)/apoln
605          tps  = SSUM(iim,tpps,1)/apols
606
607        DO ij = 1, iip1
608          ps(  ij    ) = tpn
609          ps(ij+ip1jm) = tps
610        ENDDO
611
612
613      END IF ! of IF(apdiss)
614
615c ajout debug
616c              IF( lafin ) then 
617c                abort_message = 'Simulation finished'
618c                call abort_gcm(modname,abort_message,0)
619c              ENDIF
620       
621c   ********************************************************************
622c   ********************************************************************
623c   .... fin de l'integration dynamique  et physique pour le pas itau ..
624c   ********************************************************************
625c   ********************************************************************
626
627c   preparation du pas d'integration suivant  ......
628
629      IF ( .NOT.purmats ) THEN
630c       ........................................................
631c       ..............  schema matsuno + leapfrog  ..............
632c       ........................................................
633
634            IF(forward. OR. leapf) THEN
635              itau= itau + 1
636c              iday= day_ini+itau/day_step
637c              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
638c                IF(time.GT.1.) THEN
639c                  time = time-1.
640c                  iday = iday+1
641c                ENDIF
642            ENDIF
643
644
645            IF( itau. EQ. itaufinp1 ) then 
646              if (flag_verif) then
647                write(79,*) 'ucov',ucov
648                write(80,*) 'vcov',vcov
649                write(81,*) 'teta',teta
650                write(82,*) 'ps',ps
651                write(83,*) 'q',q
652                WRITE(85,*) 'q1 = ',q(:,:,1)
653                WRITE(86,*) 'q3 = ',q(:,:,3)
654              endif
655
656              abort_message = 'Simulation finished'
657
658              call abort_gcm(modname,abort_message,0)
659            ENDIF
660c-----------------------------------------------------------------------
661c   ecriture du fichier histoire moyenne:
662c   -------------------------------------
663
664            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
665               IF(itau.EQ.itaufin) THEN
666                  iav=1
667               ELSE
668                  iav=0
669               ENDIF
670               
671               IF (ok_dynzon) THEN
672#ifdef CPP_IOIPSL
673c les traceurs ne sont pas sortis, trop lourd.
674c Peut changer eventuellement si besoin.
675                 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
676     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
677     &                 du,dudis,dutop,dufi)
678#endif
679               END IF
680               IF (ok_dyn_ave) THEN
681#ifdef CPP_IOIPSL
682                 CALL writedynav(itau,vcov,
683     &                 ucov,teta,pk,phi,q,masse,ps,phis)
684#endif
685               ENDIF
686
687            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
688
689c-----------------------------------------------------------------------
690c   ecriture de la bande histoire:
691c   ------------------------------
692
693            IF( MOD(itau,iecri).EQ.0) THEN
694             ! Ehouarn: output only during LF or Backward Matsuno
695             if (leapf.or.(.not.leapf.and.(.not.forward))) then
696              nbetat = nbetatdem
697! ADAPTATION GCM POUR CP(T)
698              call tpot2t(ijp1llm,teta,temp,pk)
699              tsurpk = cpp*temp/pk
700              CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
701              unat=0.
702              do l=1,llm
703                unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
704                vnat(:,l)=vcov(:,l)/cv(:)
705              enddo
706#ifdef CPP_IOIPSL
707              if (ok_dyn_ins) then
708!               write(lunout,*) "leapfrog: call writehist, itau=",itau
709               CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
710!               call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
711!               call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
712!              call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
713!               call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
714!               call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
715              endif ! of if (ok_dyn_ins)
716#endif
717! For some Grads outputs of fields
718              if (output_grads_dyn) then
719#include "write_grads_dyn.h"
720              endif
721             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
722            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
723
724            IF(itau.EQ.itaufin) THEN
725
726
727              if (planet_type.eq."mars") then
728! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
729                abort_message = 'dynredem1_mars A FAIRE'
730                call abort_gcm(modname,abort_message,0)
731              else
732                CALL dynredem1("restart.nc",0.0,
733     &                         vcov,ucov,teta,q,masse,ps)
734              endif ! of if (planet_type.eq."mars")
735
736              CLOSE(99)
737            ENDIF ! of IF (itau.EQ.itaufin)
738
739c-----------------------------------------------------------------------
740c   gestion de l'integration temporelle:
741c   ------------------------------------
742
743            IF( MOD(itau,iperiod).EQ.0 )    THEN
744                    GO TO 1
745            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
746
747                   IF( forward )  THEN
748c      fin du pas forward et debut du pas backward
749
750                      forward = .FALSE.
751                        leapf = .FALSE.
752                           GO TO 2
753
754                   ELSE
755c      fin du pas backward et debut du premier pas leapfrog
756
757                        leapf =  .TRUE.
758                        dt  =  2.*dtvr
759                        GO TO 2
760                   END IF ! of IF (forward)
761            ELSE
762
763c      ......   pas leapfrog  .....
764
765                 leapf = .TRUE.
766                 dt  = 2.*dtvr
767                 GO TO 2
768            END IF ! of IF (MOD(itau,iperiod).EQ.0)
769                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
770
771      ELSE ! of IF (.not.purmats)
772
773c       ........................................................
774c       ..............       schema  matsuno        ...............
775c       ........................................................
776            IF( forward )  THEN
777
778             itau =  itau + 1
779c             iday = day_ini+itau/day_step
780c             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
781c
782c                  IF(time.GT.1.) THEN
783c                   time = time-1.
784c                   iday = iday+1
785c                  ENDIF
786
787               forward =  .FALSE.
788               IF( itau. EQ. itaufinp1 ) then 
789                 abort_message = 'Simulation finished'
790                 call abort_gcm(modname,abort_message,0)
791               ENDIF
792               GO TO 2
793
794            ELSE ! of IF(forward) i.e. backward step
795
796              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
797               IF(itau.EQ.itaufin) THEN
798                  iav=1
799               ELSE
800                  iav=0
801               ENDIF
802
803               IF (ok_dynzon) THEN
804#ifdef CPP_IOIPSL
805c les traceurs ne sont pas sortis, trop lourd.
806c Peut changer eventuellement si besoin.
807                 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
808     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
809     &                 du,dudis,dutop,dufi)
810#endif
811               ENDIF
812               IF (ok_dyn_ave) THEN
813#ifdef CPP_IOIPSL
814                 CALL writedynav(itau,vcov,
815     &                 ucov,teta,pk,phi,q,masse,ps,phis)
816#endif
817               ENDIF
818
819              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
820
821              IF(MOD(itau,iecri         ).EQ.0) THEN
822c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
823                nbetat = nbetatdem
824! ADAPTATION GCM POUR CP(T)
825                call tpot2t(ijp1llm,teta,temp,pk)
826                tsurpk = cpp*temp/pk
827                CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
828                unat=0.
829                do l=1,llm
830                  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
831                  vnat(:,l)=vcov(:,l)/cv(:)
832                enddo
833#ifdef CPP_IOIPSL
834              if (ok_dyn_ins) then
835!                write(lunout,*) "leapfrog: call writehist (b)",
836!     &                        itau,iecri
837                CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
838              endif ! of if (ok_dyn_ins)
839#endif
840! For some Grads outputs
841                if (output_grads_dyn) then
842#include "write_grads_dyn.h"
843                endif
844
845              ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
846
847              IF(itau.EQ.itaufin) THEN
848                if (planet_type.eq."mars") then
849! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
850                  abort_message = 'dynredem1_mars A FAIRE'
851                  call abort_gcm(modname,abort_message,0)
852                else
853                  CALL dynredem1("restart.nc",0.0,
854     &                         vcov,ucov,teta,q,masse,ps)
855                endif ! of if (planet_type.eq."mars")
856              ENDIF ! of IF(itau.EQ.itaufin)
857
858              forward = .TRUE.
859              GO TO  1
860
861            ENDIF ! of IF (forward)
862
863      END IF ! of IF(.not.purmats)
864
865      first=.false.
866
867      STOP
868      END
Note: See TracBrowser for help on using the repository browser.