Ignore:
Timestamp:
Apr 9, 2009, 12:11:35 PM (15 years ago)
Author:
Laurent Fairhead
Message:

Réintegration dans le tronc des modifications issues de la branche LMDZ-dev
comprises entre la révision 1074 et 1145
Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes
résultats entre la trunk et la dev
LF

Location:
LMDZ4/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk

  • LMDZ4/trunk/libf/dyn3d/leapfrog.F

    r1060 r1146  
    22c
    33c
    4       SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
     4      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
    55     &                    time_0)
    66
    77
    88cIM : pour sortir les param. du modele dans un fis. netcdf 110106
    9       USE IOIPSL
     9#ifdef CPP_IOIPSL
     10      use IOIPSL
     11#endif
     12      USE infotrac
    1013
    1114      IMPLICIT NONE
    … …  
    5659#include "com_io_dyn.h"
    5760#include "iniprint.h"
    58 #include "advtrac.h"
    59 c#include "tracstoke.h"
    60 
    6161#include "academic.h"
    6262
    6363! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
    6464! #include "clesphys.h"
    65 
    66       integer nq
    6765
    6866      INTEGER         longcles
    … …  
    7674      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    7775      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    78       REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
     76      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
    7977      REAL ps(ip1jmp1)                       ! pression  au sol
    8078      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    … …  
    9795c   tendances dynamiques
    9896      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
    99       REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
     97      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
    10098
    10199c   tendances de la dissipation
    … …  
    105103c   tendances physiques
    106104      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
    107       REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
     105      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
    108106
    109107c   variables pour le fichier histoire
    … …  
    165163
    166164      character*80 dynhist_file, dynhistave_file
    167       character*20 modname
     165      character(len=20) :: modname
    168166      character*80 abort_message
    169167
    … …  
    182180      PARAMETER (testita = 9)
    183181
    184       logical , parameter :: flag_verif = .false.
     182      logical , parameter :: flag_verif = .true.
    185183     
    186184
    … …  
    190188      itaufin   = nday*day_step
    191189      itaufinp1 = itaufin +1
    192 
     190      modname="leapfrog"
     191     
    193192
    194193      itau = 0
    … …  
    220219        call guide(itau,ucov,vcov,teta,q,masse,ps)
    221220      else
    222         IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
    223      .    '6 dernieres heures'
     221        IF(prt_level>9)WRITE(lunout,*)'leapfrog: attention on ne ',
     222     .    'guide pas les 6 dernieres heures'
    224223      endif
    225224#endif
    … …  
    230229c     ENDIF
    231230c
     231
     232! Save fields obtained at previous time step as '...m1'
    232233      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
    233234      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
    … …  
    245246      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
    246247
    247       call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
     248! Ehouarn: what is this for? zqmin & zqmax are not used anyway ...
     249!      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
    248250
    249251   2  CONTINUE
    … …  
    305307
    306308
    307          ENDIF
    308 c
    309       ENDIF
     309         ENDIF ! of IF (offline)
     310c
     311      ENDIF ! of IF( forward. OR . leapf )
    310312
    311313
    … …  
    353355c   -----------------------------------------------------
    354356
    355 #ifdef CPP_PHYS
    356357c+jld
    357358
    358359c  Diagnostique de conservation de l'énergie : initialisation
    359       IF (ip_ebil_dyn.ge.1 ) THEN
     360         IF (ip_ebil_dyn.ge.1 ) THEN
    360361          ztit='bil dyn'
    361           CALL diagedyn(ztit,2,1,1,dtphys
    362      e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
    363       ENDIF
     362! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
     363           IF (planet_type.eq."earth") THEN
     364            CALL diagedyn(ztit,2,1,1,dtphys
     365     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
     366           ENDIF
     367         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
    364368c-jld
     369#ifdef CPP_IOIPSL
    365370cIM : pour sortir les param. du modele dans un fis. netcdf 110106
    366       IF (first) THEN
    367        first=.false.
     371         IF (first) THEN
     372          first=.false.
    368373#include "ini_paramLMDZ_dyn.h"
    369       ENDIF
     374         ENDIF
    370375c
    371376#include "write_paramLMDZ_dyn.h"
    372377c
    373 
    374         CALL calfis( nq, lafin ,rdayvrai,time  ,
     378#endif
     379! #endif of #ifdef CPP_IOIPSL
     380         CALL calfis( lafin ,rdayvrai,time  ,
    375381     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
    376382     $               du,dv,dteta,dq,
    … …  
    378384     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
    379385
    380        IF (ok_strato) THEN
    381          CALL top_bound( vcov,ucov,teta, dufi,dvfi,dtetafi)
    382        ENDIF
     386         IF (ok_strato) THEN
     387           CALL top_bound( vcov,ucov,teta, dufi,dvfi,dtetafi)
     388         ENDIF
    383389       
    384390c      ajout des tendances physiques:
    385391c      ------------------------------
    386           CALL addfi( nqmx, dtphys, leapf, forward   ,
     392          CALL addfi( dtphys, leapf, forward   ,
    387393     $                  ucov, vcov, teta , q   ,ps ,
    388394     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    389395c
    390396c  Diagnostique de conservation de l'énergie : difference
    391       IF (ip_ebil_dyn.ge.1 ) THEN
     397         IF (ip_ebil_dyn.ge.1 ) THEN
    392398          ztit='bil phys'
    393           CALL diagedyn(ztit,2,1,1,dtphys
    394      e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
    395       ENDIF
    396 #endif
     399          IF (planet_type.eq."earth") THEN
     400           CALL diagedyn(ztit,2,1,1,dtphys
     401     &     , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
     402          ENDIF
     403         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
     404
    397405       ENDIF ! of IF( apphys )
    398406
    399       IF(iflag_phys.EQ.2) THEN ! "Newtonian physics" case
     407      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
    400408c   Calcul academique de la physique = Rappel Newtonien + friction
    401409c   --------------------------------------------------------------
    … …  
    475483
    476484
    477       END IF
     485      END IF ! of IF(apdiss)
    478486
    479487c ajout debug
    … …  
    509517            IF( itau. EQ. itaufinp1 ) then 
    510518              if (flag_verif) then
    511                 write(79,*) 'ucov',ucov
    512                 write(80,*) 'vcov',vcov
    513                 write(81,*) 'teta',teta
    514                 write(82,*) 'ps',ps
    515                 write(83,*) 'q',q
     519                write(80,*) 'ucov',ucov
     520                write(81,*) 'vcov',vcov
     521                write(82,*) 'teta',teta
     522                write(83,*) 'ps',ps
     523                write(84,*) 'q',q
    516524                WRITE(85,*) 'q1 = ',q(:,:,1)
    517525                WRITE(86,*) 'q3 = ',q(:,:,3)
     526                write(90) ucov
     527                write(91) vcov
     528                write(92) teta
     529                write(93) ps
     530                write(94) q
    518531              endif
    519532
    … …  
    532545                  iav=0
    533546               ENDIF
     547               
     548               IF (ok_dynzon) THEN
    534549#ifdef CPP_IOIPSL
    535               CALL writedynav(histaveid, nqmx, itau,vcov ,
    536      ,                          ucov,teta,pk,phi,q,masse,ps,phis)
    537                call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
    538      ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
    539 #endif
     550                  CALL writedynav(histaveid, itau,vcov ,
     551     ,                 ucov,teta,pk,phi,q,masse,ps,phis)
     552                  CALL bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
     553     ,                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
     554#endif
     555               END IF
    540556
    541557            ENDIF
    … …  
    548564c           IF( MOD(itau,iecri*day_step).EQ.0) THEN
    549565
    550                nbetat = nbetatdem
    551        CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
    552         unat=0.
    553         do l=1,llm
    554            unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
    555            vnat(:,l)=vcov(:,l)/cv(:)
    556         enddo
     566              nbetat = nbetatdem
     567              CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     568              unat=0.
     569              do l=1,llm
     570                unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
     571                vnat(:,l)=vcov(:,l)/cv(:)
     572              enddo
    557573#ifdef CPP_IOIPSL
    558 c        CALL writehist(histid,histvid, nqmx,itau,vcov,
    559 c     s                       ucov,teta,phi,q,masse,ps,phis)
    560 #else
     574c             CALL writehist(histid,histvid,itau,vcov,
     575c     &                      ucov,teta,phi,q,masse,ps,phis)
     576#endif
     577! For some Grads outputs of fields
     578             if (output_grads_dyn) then
    561579#include "write_grads_dyn.h"
    562 #endif
    563 
    564 
    565             ENDIF
     580             endif
     581
     582            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
    566583
    567584            IF(itau.EQ.itaufin) THEN
    568585
    569586
    570 #ifdef CPP_IOIPSL
    571        CALL dynredem1("restart.nc",0.0,
    572      ,                     vcov,ucov,teta,q,nqmx,masse,ps)
    573 #endif
     587              if (planet_type.eq."earth") then
     588#ifdef CPP_EARTH
     589! Write an Earth-format restart file
     590                CALL dynredem1("restart.nc",0.0,
     591     &                         vcov,ucov,teta,q,masse,ps)
     592#endif
     593              endif ! of if (planet_type.eq."earth")
    574594
    575595              CLOSE(99)
    576             ENDIF
     596            ENDIF ! of IF (itau.EQ.itaufin)
    577597
    578598c-----------------------------------------------------------------------
    … …  
    596616                        leapf =  .TRUE.
    597617                        dt  =  2.*dtvr
    598                         GO TO 2
    599                    END IF
     618                        GO TO 2 
     619                   END IF ! of IF (forward)
    600620            ELSE
    601621
    … …  
    605625                 dt  = 2.*dtvr
    606626                 GO TO 2
    607             END IF
    608 
    609       ELSE
     627            END IF ! of IF (MOD(itau,iperiod).EQ.0)
     628                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
     629
     630      ELSE ! of IF (.not.purmats)
    610631
    611632c       ........................................................
    … …  
    630651               GO TO 2
    631652
    632             ELSE
    633 
    634             IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
     653            ELSE ! of IF(forward)
     654
     655              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
    635656               IF(itau.EQ.itaufin) THEN
    636657                  iav=1
    … …  
    638659                  iav=0
    639660               ENDIF
     661
     662               IF (ok_dynzon) THEN
    640663#ifdef CPP_IOIPSL
    641               CALL writedynav(histaveid, nqmx, itau,vcov ,
    642      ,                          ucov,teta,pk,phi,q,masse,ps,phis)
    643                call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
    644      ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
    645 #endif
    646 
    647             ENDIF
    648 
    649                IF(MOD(itau,iecri         ).EQ.0) THEN
     664                  CALL writedynav(histaveid, itau,vcov ,
     665     ,                 ucov,teta,pk,phi,q,masse,ps,phis)
     666                  CALL bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
     667     ,                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
     668#endif
     669               END IF
     670
     671              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
     672
     673              IF(MOD(itau,iecri         ).EQ.0) THEN
    650674c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
    651                   nbetat = nbetatdem
    652        CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
    653         unat=0.
    654         do l=1,llm
    655            unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
    656            vnat(:,l)=vcov(:,l)/cv(:)
    657         enddo
     675                nbetat = nbetatdem
     676                CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     677                unat=0.
     678                do l=1,llm
     679                  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
     680                  vnat(:,l)=vcov(:,l)/cv(:)
     681                enddo
    658682#ifdef CPP_IOIPSL
    659 c       CALL writehist( histid, histvid, nqmx, itau,vcov ,
    660 c    ,                           ucov,teta,phi,q,masse,ps,phis)
    661 #else
     683c               CALL writehist( histid, histvid, itau,vcov ,
     684c    &                           ucov,teta,phi,q,masse,ps,phis)
     685#endif
     686! For some Grads outputs
     687                if (output_grads_dyn) then
    662688#include "write_grads_dyn.h"
    663 #endif
    664 
    665 
    666                ENDIF
    667 
    668 #ifdef CPP_IOIPSL
    669                  IF(itau.EQ.itaufin)
    670      . CALL dynredem1("restart.nc",0.0,
    671      .                     vcov,ucov,teta,q,nqmx,masse,ps)
    672 #endif
    673 
    674                  forward = .TRUE.
    675                  GO TO  1
    676 
    677             ENDIF
    678 
    679       END IF
     689                endif
     690
     691              ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
     692
     693              IF(itau.EQ.itaufin) THEN
     694                if (planet_type.eq."earth") then
     695#ifdef CPP_EARTH
     696                  CALL dynredem1("restart.nc",0.0,
     697     &                           vcov,ucov,teta,q,masse,ps)
     698#endif
     699                endif ! of if (planet_type.eq."earth")
     700              ENDIF ! of IF(itau.EQ.itaufin)
     701
     702              forward = .TRUE.
     703              GO TO  1
     704
     705            ENDIF ! of IF (forward)
     706
     707      END IF ! of IF(.not.purmats)
    680708
    681709      STOP
Note: See TracChangeset for help on using the changeset viewer.