Ignore:
Timestamp:
Dec 9, 2024, 3:39:31 PM (16 months ago)
Author:
tbertrand
Message:

LMDZ.PLUTO:
Update and Validation of the Volatile Transport Model (VTM, or "nogcm")
Merging missing elements with Pluto.old
TB

Location:
trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/leapfrog_nogcm.F

    r3422 r3539  
    3838      IMPLICIT NONE
    3939
    40 c      ......   Version  du 10/01/98    ..........
    41 
    42 c             avec  coordonnees  verticales hybrides
    43 c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
     40c      ......   Version  du 30/11/24    ..........
    4441
    4542c=======================================================================
    4643c
    47 c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
     44c   Auteur:  T. Bertrand, F. Forget, A. Falco
    4845c   -------
    4946c
    … …  
    5148c   ------
    5249c
    53 c   GCM LMD nouvelle grille
     50c   GCM LMD sans dynamique, pour modele saisonnier de surface
    5451c
    5552c=======================================================================
    56 c
    57 c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
    58 c      et possibilite d'appeler une fonction f(y)  a derivee tangente
    59 c      hyperbolique a la  place de la fonction a derivee sinusoidale.
    60 
    61 c  ... Possibilite de choisir le shema pour l'advection de
    62 c        q  , en modifiant iadv dans traceur.def  (10/02) .
    63 c
    64 c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
    65 c      Pour Van-Leer iadv=10
    6653c
    6754c-----------------------------------------------------------------------
    … …  
    7663#include "iniprint.h"
    7764#include "academic.h"
    78 
    79 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
    80 ! #include "clesphys.h"
    8165
    8266      REAL,INTENT(IN) :: time_0 ! not used
    … …  
    10387      REAL tsurpk(ip1jmp1,llm)               ! cpp*T/pk 
    10488
    105 !      real zqmin,zqmax
    106 
    107 
    108 !       ED18 nogcm
     89!     nogcm
    10990      REAL tau_ps
    11091      REAL tau_x
    … …  
    11798      REAL tetamean
    11899      real globaverage2d
    119 !      LOGICAL firstcall_globaverage2d
    120 
    121100
    122101c variables dynamiques intermediaire pour le transport
    … …  
    132111      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot)
    133112
    134 c   tendances de la dissipation en */s
    135       REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
    136       REAL dtetadis(ip1jmp1,llm)
    137 
    138 c   tendances de la couche superieure */s
    139 c      REAL dvtop(ip1jm,llm)
    140       REAL dutop(ip1jmp1,llm)
    141 c      REAL dtetatop(ip1jmp1,llm)
    142 c      REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
    143 
    144 c   TITAN : tendances due au forces de marees */s
    145       REAL dvtidal(ip1jm,llm),dutidal(ip1jmp1,llm)
    146 
    147113c   tendances physiques */s
    148114      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
    … …  
    159125
    160126      REAL  SSUM
    161 !     REAL finvmaold(ip1jmp1,llm)
    162127
    163128cym      LOGICAL  lafin
    164129      LOGICAL :: lafin=.false.
    165130      INTEGER ij,iq,l
    166 !      INTEGER ik
    167 
    168 !      real time_step, t_wrt, t_ops
    169131
    170132      REAL rdaym_ini
    … …  
    179141      save first
    180142      data first/.true./
    181 !      real dt_cum
    182 !      character*10 infile
    183 !      integer zan, tau0, thoriid
    184 !      integer nid_ctesGCM
    185 !      save nid_ctesGCM
    186 !      real degres
    187 !      real rlong(iip1), rlatg(jjp1)
    188 !      real zx_tmp_2d(iip1,jjp1)
    189 !      integer ndex2d(iip1*jjp1)
    190143      logical ok_sync
    191144      parameter (ok_sync = .true.)
    … …  
    214167      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
    215168      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
    216 !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
    217169      CHARACTER*15 ztit
    218 !IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
    219 !IM   SAVE      ip_ebil_dyn
    220 !IM   DATA      ip_ebil_dyn/0/
    221 c-jld
    222 
    223 !      integer :: itau_w ! for write_paramLMDZ_dyn.h
    224 
    225 !      character*80 dynhist_file, dynhistave_file
     170
    226171      character(len=*),parameter :: modname="leapfrog"
    227172      character*80 abort_message
    … …  
    271216
    272217c INITIALISATIONS
    273         dudis(:,:)   =0.
    274         dvdis(:,:)   =0.
    275         dtetadis(:,:)=0.
    276         dutop(:,:)   =0.
    277 c        dvtop(:,:)   =0.
    278 c        dtetatop(:,:)=0.
    279 c        dqtop(:,:,:) =0.
    280 c        dptop(:)     =0.
    281218        dufi(:,:)   =0.
    282219        dvfi(:,:)   =0.
    … …  
    330267      write(*,*) "Fin Test PK"
    331268c     stop
     269
    332270c------------------
    333271c     Preparing mixing of pressure and tracers in nogcm
    … …  
    337275c     It is checked that globalaverage2d(Pi)=pmean
    338276      DO ij=1,ip1jmp1
    339              kpd(ij) = exp(-phis(ij)/(r*200.))
     277             kpd(ij) = exp(-phis(ij)/(r*38.))
    340278      ENDDO
    341279      p00d=globaverage2d(kpd) ! mean pres at ref level
    342       tau_ps = 1. ! constante de rappel for pressure  (s)
    343       tau_n2 = 1 ! constante de rappel for mix ratio qn2 (s)
    344       tau_def = 1.E7 ! default constante de rappel for mix ratio qX (s)
     280      tau_ps = tau_n2 ! constante de rappel for pressure  (s)
     281      tau_def = 1. ! default constante de rappel for mix ratio qX (s)
    345282      tau_teta = 1.E7 !constante de rappel for potentiel temperature
    346283
    … …  
    368305c   -----
    369306
    370       jD_cur = jD_ref + day_ini - day_ref +                             &
     307      jD_cur = jD_ref + day_ini - day_ref +                             
    371308     &          (itau+1)/day_step
    372       jH_cur = jH_ref + start_time +                                    &
     309      jH_cur = jH_ref + start_time +                                   
    373310     &          mod(itau+1,day_step)/float(day_step)
    374311      jD_cur = jD_cur + int(jH_cur)
    375312      jH_cur = jH_cur - int(jH_cur)
    376 
    377 c
    378313
    379314! Save fields obtained at previous time step as '...m1'
    … …  
    407342c   gestion des appels de la physique et des dissipations:
    408343c   ------------------------------------------------------
    409 c
    410 c   ...    P.Le Van  ( 6/02/95 )  ....
    411 
    412 ! ED18: suppression des mentions de la variable apdiss dans le cas
    413 ! 'nogcm'
    414344
    415345      apphys = .FALSE.
    416346      statcl = .FALSE.
    417347      conser = .FALSE.
    418      
    419348
    420349      IF( purmats ) THEN
    421350      ! Purely Matsuno time stepping
    422351         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
    423    
    424352   
    425353         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
    426354     s          .and. physics        ) apphys = .TRUE.
    427355      ELSE
    428       ! Leapfrog/Matsuno time stepping
    429         IF( MOD(itau   ,iconser) .EQ. 0  ) conser = .TRUE.
    430  
    431  
    432         IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
     356         ! Leapfrog/Matsuno time stepping
     357         IF( MOD(itau   ,iconser) .EQ. 0  ) conser = .TRUE.
     358         IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
    433359      END IF
    434 
    435 ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
    436 !          supress dissipation step
    437 
    438 
    439360
    440361c-----------------------------------------------------------------------
    … …  
    442363c   --------------------------------
    443364
    444 ! ED18: suppression de l'onglet pour le cas nogcm
    445 
    446 
    447          dv(:,:) = 0.
    448          du(:,:) = 0.
    449          dteta(:,:) = 0.
    450          dq(:,:,:) = 0.
    451          
    452 
    453 
    454 c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
    455 c
     365      dv(:,:) = 0.
     366      du(:,:) = 0.
     367      dteta(:,:) = 0.
     368      dq(:,:,:) = 0.
     369
    456370c-----------------------------------------------------------------------
    457371c   calcul des tendances physiques:
    458372c   -------------------------------
    459 c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
    460373c
    461374       IF( purmats )  THEN
    … …  
    464377          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
    465378       ENDIF
    466 c
    467 c
     379
    468380       IF( apphys )  THEN
    469 c
    470 c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
    471 c
    472381
    473382         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
    … …  
    481390         CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
    482391
    483            jD_cur = jD_ref + day_ini - day_ref +                        &
     392         jD_cur = jD_ref + day_ini - day_ref +                       
    484393     &          (itau+1)/day_step
    485394
    486            ! AS: we make jD_cur to be pday
    487            jD_cur = int(day_ini + itau/day_step)
    488 
    489 !           print*,'itau =',itau
    490 !           print*,'day_step =',day_step
    491 !           print*,'itau/day_step =',itau/day_step
    492 
    493 
    494            jH_cur = jH_ref + start_time +                               &
     395         ! AS: we make jD_cur to be pday
     396         jD_cur = int(day_ini + itau/day_step)
     397
     398         jH_cur = jH_ref + start_time +                               
    495399     &          mod(itau+1,day_step)/float(day_step)
    496            jH_cur = jH_ref + start_time +                               &
     400         jH_cur = jH_ref + start_time +                               
    497401     &          mod(itau,day_step)/float(day_step)
    498            jD_cur = jD_cur + int(jH_cur)
    499            jH_cur = jH_cur - int(jH_cur)
    500 
    501            
    502 
    503 !         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
    504 !         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
    505 !         write(lunout,*)'current date = ',an, mois, jour, secondes
    506 
    507 c rajout debug
    508 c       lafin = .true.
    509 
     402         jD_cur = jD_cur + int(jH_cur)
     403         jH_cur = jH_cur - int(jH_cur)
    510404
    511405c   Interface avec les routines de phylmd (phymars ... )
    512406c   -----------------------------------------------------
    513 
    514 c+jld
    515 
    516 c  Diagnostique de conservation de l'Energie : initialisation
    517          IF (ip_ebil_dyn.ge.1 ) THEN
    518           ztit='bil dyn'
    519 ! Ehouarn: be careful, diagedyn is Earth-specific!
    520            IF (planet_type.eq."earth") THEN
    521             CALL diagedyn(ztit,2,1,1,dtphys
    522      &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
    523            ENDIF
    524          ENDIF ! of IF (ip_ebil_dyn.ge.1 )
    525 c-jld
    526 #ifdef CPP_IOIPSL
    527 cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
    528 cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
    529 c        IF (first) THEN
    530 c         first=.false.
    531 c#include "ini_paramLMDZ_dyn.h"
    532 c        ENDIF
    533 c
    534 c#include "write_paramLMDZ_dyn.h"
    535 c
    536 #endif
    537 ! #endif of #ifdef CPP_IOIPSL
    538 
    539 c          call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
    540 !         print*,'---------LEAPFROG---------------'
    541 !         print*,''
    542 !         print*,'> AVANT CALFIS'
    543 !         print*,''
    544 !         print*,'teta(3049,:) = ',teta(3049,:)
    545 !         print*,''
    546 !         print*,'dteta(3049,:) = ',dteta(3049,:)
    547 !         print*,''
    548 !         print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
    549 !         print*,''
    550407
    551408         CALL calfis( lafin , jD_cur, jH_cur,
    … …  
    555412     $               dufi,dvfi,dtetafi,dqfi,dpfi  )
    556413
    557 c          call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
    558 c          call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
    559 c          call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
    560 !         print*,'> APRES CALFIS (AVANT ADDFI)'
    561 !         print*,''
    562 !         print*,'teta(3049,:) = ',teta(3049,:)
    563 !         print*,''
    564 !         print*,'dteta(3049,:) = ',dteta(3049,:)
    565 !         print*,''
    566 !         print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
    567 !         print*,''
    568 
    569 
    570414c      ajout des tendances physiques
    571415c      ------------------------------
    572416
    573           CALL addfi( dtphys, leapf, forward   ,
     417         CALL addfi( dtphys, leapf, forward   ,
    574418     $                  ucov, vcov, teta , q   ,ps ,
    575419     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    576420
    577 !         print*,'> APRES ADDFI'
    578 !         print*,''
    579 !         print*,'teta(3049,:) = ',teta(3049,:)
    580 !         print*,''
    581 !         print*,'dteta(3049,:) = ',dteta(3049,:)
    582 !         print*,''
    583 !         print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
    584 !         print*,''
    585 
    586           CALL pression (ip1jmp1,ap,bp,ps,p)
    587           CALL massdair(p,masse)
     421         CALL pression (ip1jmp1,ap,bp,ps,p)
     422         CALL massdair(p,masse)
    588423   
    589          ! Mass of tracers in each mess (kg) before Horizontal mixing of pressure
    590 c        -------------------------------------------------------------------
     424         ! 1) Mass of tracers in each mess (kg) before Horizontal mixing of pressure
    591425         
    592426         DO l=1, llm
    … …  
    596430         ENDDO
    597431
    598 c        Horizontal mixing of pressure
    599 c        ------------------------------
     432         ! 2) Horizontal mixing of pressure
    600433         ! Rappel newtonien vers psmean
    601            psmean= globaverage2d(ps)  ! mean pressure
    602            p0=psmean/p00d
    603            DO ij=1,ip1jmp1
    604                 oldps(ij)=ps(ij)
    605                 ps(ij)=ps(ij) +(p0*kpd(ij)
     434         psmean= globaverage2d(ps)  ! mean pressure
     435         p0=psmean/p00d
     436         DO ij=1,ip1jmp1
     437             oldps(ij)=ps(ij)
     438             ps(ij)=ps(ij) +(p0*kpd(ij)
    606439     &                 -ps(ij))*(1-exp(-dtphys/tau_ps))
    607            ENDDO
    608 
    609            write(72,*) 'psmean = ',psmean  ! mean pressure redistributed
    610 
    611          ! security check: test pressure negative
    612            DO ij=1,ip1jmp1
     440         ENDDO
     441
     442         write(72,*) 'psmean = ',psmean  ! mean pressure redistributed
     443
     444         ! 3) Security check: test pressure negative
     445         DO ij=1,ip1jmp1
    613446             IF (ps(ij).le.0.) then
    614447                 !write(*,*) 'Negative pressure :'
    … …  
    617450                 !write(*,*) 'tau_ps = ',tau_ps
    618451                 !STOP
    619                  ps(ij)=0.0000001*kpd(ij)/p00d
     452                 ps(ij)=1.e-6*kpd(ij)/p00d
    620453             ENDIF
    621            ENDDO
    622 !***********************
    623 !          Correction on teta due to surface pressure changes
    624            DO l = 1,llm
    625             DO ij = 1,ip1jmp1
    626               teta(ij,l)= teta(ij,l)*(ps(ij)/oldps(ij))**kappa
    627             ENDDO
    628            ENDDO
    629 !***********************
    630 
    631 
    632 
    633 !        ! update pressure and update p and pk
    634 !         DO ij=1,ip1jmp1
    635 !            ps(ij) = ps(ij) + dpfi(ij)*dtphys
     454         ENDDO
     455
     456!         ! 4) Correction on teta due to surface pressure changes
     457!         DO l = 1,llm
     458!            DO ij = 1,ip1jmp1
     459!              teta(ij,l)= teta(ij,l)*(ps(ij)/oldps(ij))**kappa
     460!            ENDDO
    636461!         ENDDO
    637           CALL pression (ip1jmp1,ap,bp,ps,p)
    638           CALL massdair(p,masse)
    639           if (pressure_exner) then
     462
     463         ! 5) Update pressure p and pk
     464         CALL pression (ip1jmp1,ap,bp,ps,p)
     465         CALL massdair(p,masse)
     466         if (pressure_exner) then
    640467            CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
    641           else
     468         else
    642469            CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    643           endif
    644 
    645          ! Update tracers after Horizontal mixing of pressure to ! conserve  tracer mass
    646 c        -------------------------------------------------------------------
     470         endif
     471
     472         ! 6) Update tracers after Horizontal mixing of pressure to ! conserve  tracer mass
    647473         DO l=1, llm
    648474           DO ij=1,ip1jmp1
    … …  
    650476           ENDDO
    651477         ENDDO
    652          
    653 
    654 
    655 c        Horizontal mixing of pressure
    656 c        ------------------------------
    657          ! Rappel newtonien vers psmean
    658            psmean= globaverage2d(ps)  ! mean pressure
    659 !        ! increment q_n2  with physical tendancy
    660 !          IF (igcm_n2.ne.0) then
    661 !            DO l=1, llm
    662 !               DO ij=1,ip1jmp1
    663 !                q(ij,l,igcm_n2)=q(ij,l,igcm_n2)+
    664 !    &                    dqfi(ij,l,igcm_n2)*dtphys
    665 !               ENDDO
    666 !            ENDDO
    667 !          ENDIF
    668 
    669 c          Mixing N2 vertically ! not used for pluto ?
    670 c          --------------------------
    671 !            if (igcm_n2.ne.0) then
    672 !             DO ij=1,ip1jmp1
    673 !                qmean_x_vert=0.
    674 !                DO l=1, llm
    675 !                  qmean_x_vert= qmean_x_vert
    676 !      &          + q(ij,l,igcm_n2)*( p(ij,l) - p(ij,l+1))
    677 !                END DO
    678 !                qmean_x_vert= qmean_x_vert/ps(ij)
    679 !                DO l=1, llm
    680 !                  q(ij,l,igcm_n2)= qmean_x_vert
    681 !                END DO
    682 !             END DO
    683 !            end if
    684 
    685 
    686 
    687 c        Horizontal mixing of pressure
    688 c        ------------------------------
    689          ! Rappel newtonien vers psmean
    690            psmean= globaverage2d(ps)  ! mean pressure
    691 
    692 c        Horizontal mixing tracers, and temperature (replace dynamics in nogcm)
    693 c        --------------------------------------------------------------------------------- 
     478
     479         ! 7) Horizontal mixing tracers, and temperature (replace dynamics in nogcm)
     480         psmean= globaverage2d(ps)  ! mean pressure
    694481
    695482         ! Simulate redistribution by dynamics for qX
    … …  
    703490              pqxmean=globaverage2d(pqx)
    704491
    705          !    Rappel newtonien vers qx_mean
     492              ! Rappel newtonien vers qx_mean
    706493              qmean_x= pqxmean / psmean
    707494             
    … …  
    728515              END DO
    729516
    730 !             TEMPORAIRE (ED)
    731 !            PRINT*,'psmean = ',psmean
    732 !            PRINT*,'qmean_x = ',qmean_x
    733 !            PRINT*,'pqxmean = ',pqxmean
     517            ! Diagnostics
     518            ! PRINT*,'psmean = ',psmean
     519            ! PRINT*,'qmean_x = ',qmean_x
     520            ! PRINT*,'pqxmean = ',pqxmean
    734521            ! PRINT*,' q(50,1) = ',iq,q(50,1,iq)
    735522            ! PRINT*,' q(50,2) = ',iq,q(50,2,iq)
    … …  
    737524
    738525           endif ! igcm_n2.ne.0
    739       enddo
    740 
    741 
    742 !       *****************************************s***************
    743 
    744 c        Horizontal mixing of pressure
    745 c        ------------------------------
    746          ! Rappel newtonien vers psmean
    747            psmean= globaverage2d(ps)  ! mean pressure
    748 
    749 
    750 c          Mixing Temperature horizontally
    751 c          -------------------------------
    752 !          initialize variables that will be averaged
    753            DO l=1,llm
    754              DO ij=1,ip1jmp1
    755                dp(ij,l) = p(ij,l) - p(ij,l+1)
    756                tetadp(ij,l) = teta(ij,l)*dp(ij,l)
    757              ENDDO
    758            ENDDO
    759 
    760            DO l=1,llm
    761              tetadpmean = globaverage2d(tetadp(:,l))
    762              dpmean = globaverage2d(dp(:,l))
    763              tetamean = tetadpmean / dpmean
    764              DO ij=1,ip1jmp1
    765                teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
    766      &                      (1 - exp(-dtphys/tau_teta))
    767              ENDDO
    768            ENDDO
     526         ENDDO
     527
     528         ! 8) Mixing Temperature horizontally
     529         ! initialize variables that will be averaged
     530!         DO l=1,llm
     531!             DO ij=1,ip1jmp1
     532!               dp(ij,l) = p(ij,l) - p(ij,l+1)
     533!               tetadp(ij,l) = teta(ij,l)*dp(ij,l)
     534!             ENDDO
     535!         ENDDO
     536
     537!         DO l=1,llm
     538!             tetadpmean = globaverage2d(tetadp(:,l))
     539!             dpmean = globaverage2d(dp(:,l))
     540!             tetamean = tetadpmean / dpmean
     541!             DO ij=1,ip1jmp1
     542!               teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
     543!     &                      (1 - exp(-dtphys/tau_teta))
     544!             ENDDO
     545!         ENDDO
    769546           
    770 
    771547       ENDIF ! of IF( apphys )
    772548
    … …  
    777553c   ********************************************************************
    778554c   ********************************************************************
    779 
    780555
    781556      IF ( .NOT.purmats ) THEN
    … …  
    793568c                ENDIF
    794569            ENDIF
    795 
    796570
    797571            IF( itau. EQ. itaufinp1 ) then 
    … …  
    810584              call abort_gcm(modname,abort_message,0)
    811585            ENDIF
     586
    812587c-----------------------------------------------------------------------
    813588c   ecriture du fichier histoire moyenne:
    … …  
    824599               CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    825600
    826                IF (ok_dynzon) THEN
    827 #ifdef CPP_IOIPSL
    828 c les traceurs ne sont pas sortis, trop lourd.
    829 c Peut changer eventuellement si besoin.
    830                  CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
    831      &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
    832      &                 du,dudis,dutop,dufi)
    833 #endif
    834                END IF
    835601               IF (ok_dyn_ave) THEN
    836602#ifdef CPP_IOIPSL
    … …  
    963729               CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    964730
    965                IF (ok_dynzon) THEN
    966 #ifdef CPP_IOIPSL
    967 c les traceurs ne sont pas sortis, trop lourd.
    968 c Peut changer eventuellement si besoin.
    969                  CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
    970      &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
    971      &                 du,dudis,dutop,dufi)
    972 #endif
    973                ENDIF
    974731               IF (ok_dyn_ave) THEN
    975732#ifdef CPP_IOIPSL

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F

    r3438 r3539  
    100100      REAL tsurf(ngridmx)        ! surface temperature
    101101      REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
     102      REAL,ALLOCATABLE :: inertiedat_simple(:,:) ! thermal inertia tmp for dynamico
    102103      REAL emis(ngridmx)        ! surface emissivity
    103104      real emisread             ! added by RW
    … …  
    235236      call getin_p("nsoilmx",nsoilmx)
    236237     
     238      allocate(inertiedat_simple(ngridmx,nsoilmx))
    237239      allocate(tsoil(ngridmx,nsoilmx))
    238240      allocate(field_inputs(ngridmx,nsoilmx))
    … …  
    379381        CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
    380382     .        nqtot,day_ini,time,
    381      .        tsurf,tsoil,emis,q2,qsurf)    !) ! temporary modif by RDW
     383     .        tsurf,tsoil,emis,q2,qsurf,inertiedat_simple)    !) ! temporary modif by RDW
    382384
    383385        ! copy albedo and soil thermal inertia on (local) physics grid
    … …  
    40134015     &              llm,
    40144016     &              nqtot,dtphys,real(day_ini),0.0,
    4015      &              cell_area,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
     4017     &              cell_area,albfi,zmea,zstd,zsig,zgam,zthe)
    40164018      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
    40174019     &                dtphys,real(day_ini),
    4018      &                tsurf,tsoil,emis,q2,qsurf)
     4020     &                tsurf,tsoil,ithfi,emis,q2,qsurf)
    40194021!     &                cloudfrac,totalfrac,hice,
    40204022!     &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)

  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/nogcm.F90

    r3412 r3539  
    11!
    2 ! nogcm for Pluto, based on mars nogcm, 07/2024
    3 ! Author: A. Falco
     2! nogcm for Pluto, based on mars nogcm, 12/2024
     3! Author: A. Falco, T. Bertrand
    44!
    55!
    … …  
    4747  IMPLICIT NONE
    4848
    49   !      ......   Version  du 10/01/98    ..........
    50 
    51   !             avec  coordonnees  verticales hybrides
    52   !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
    53 
    54   !=======================================================================
    55   !
    56   !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
    57   !   -------
    58   !
    59   !   Objet:
    60   !   ------
    61   !
    62   !   GCM LMD nouvelle grille
    63   !
    64   !=======================================================================
    65   !
    66   !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
    67   !      et possibilite d'appeler une fonction f(y)  a derivee tangente
    68   !      hyperbolique a la  place de la fonction a derivee sinusoidale.
    69   !  ... Possibilite de choisir le schema pour l'advection de
    70   !        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
    71   !
    72   !      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
    73   !      Pour Van-Leer iadv=10
    74   !
    7549  !-----------------------------------------------------------------------
    7650  !   Declarations:
    … …  
    11488  LOGICAL first
    11589
    116   !      LOGICAL call_iniphys
    117   !      data call_iniphys/.true./
    118 
    119   !+jld variables test conservation energie
    120   !      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
    121   !     Tendance de la temp. potentiel d (theta)/ d t due a la
    122   !     tansformation d'energie cinetique en energie thermique
    123   !     cree par la dissipation
    12490  REAL dhecdt(ip1jmp1,llm)
    125   !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
    126   !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
    12791  CHARACTER (len=15) :: ztit
    128   !-jld
    12992
    13093
    … …  
    171134
    172135
    173 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    174 ! FH 2008/05/02
    175 ! A nettoyer. On ne veut qu'une ou deux routines d'interface
    176 ! dynamique -> physique pour l'initialisation
    177 !#ifdef CPP_PHYS
    178 !  CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
    179 !!      call initcomgeomphy ! now done in iniphysiq
    180 !#endif
    181 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    182136!
    183137! Initialisations pour Cp(T) Venus
    … …  
    194148    call ioconf_calendar('360d')
    195149    write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
    196   else if (calend == 'earth_365d') then
    197     call ioconf_calendar('noleap')
    198     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
    199   else if (calend == 'earth_366d') then
    200     call ioconf_calendar('gregorian')
    201     write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
    202   else if (calend == 'titan') then
    203 !        call ioconf_calendar('titan')
    204     write(lunout,*)'CALENDRIER CHOISI: Titan'
    205     abort_message = 'A FAIRE...'
    206     call abort_gcm(modname,abort_message,1)
    207   else if (calend == 'venus') then
    208 !        call ioconf_calendar('venus')
    209     write(lunout,*)'CALENDRIER CHOISI: Venus'
    210     abort_message = 'A FAIRE...'
    211     call abort_gcm(modname,abort_message,1)
    212   else
    213     abort_message = 'Mauvais choix de calendrier'
    214     call abort_gcm(modname,abort_message,1)
    215   endif
    216 #endif
    217 !-----------------------------------------------------------------------
    218   !
    219   !
     150  endif
     151#endif
     152 
    220153  !------------------------------------
    221154  !   Initialisation partie parallele
    222155  !------------------------------------
    223156
    224   !
    225   !
    226157  !-----------------------------------------------------------------------
    227158  !   Initialisation des traceurs
    … …  
    257188     endif
    258189 
    259      !       write(73,*) 'ucov',ucov
    260      !       write(74,*) 'vcov',vcov
    261      !       write(75,*) 'teta',teta
    262      !       write(76,*) 'ps',ps
    263      !       write(77,*) 'q',q
    264 
    265190  endif ! of if (read_start)
    266191
    … …  
    268193  ! le cas echeant, creation d un etat initial
    269194  IF (prt_level > 9) WRITE(lunout,*) &
    270        'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
     195       'NOGCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
    271196  if (.not.read_start) then
    272197     CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    … …  
    330255  ENDIF
    331256
    332 !      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
    333 !        write(lunout,*)
    334 !     .  'GCM: Attention les dates initiales lues dans le fichier'
    335 !        write(lunout,*)
    336 !     .  ' restart ne correspondent pas a celles lues dans '
    337 !        write(lunout,*)' gcm.def'
    338 !        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
    339 !        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
    340 !        if (raz_date .ne. 1) then
    341 !          write(lunout,*)
    342 !     .    'GCM: On garde les dates du fichier restart'
    343 !        else
    344 !          annee_ref = anneeref
    345 !          day_ref = dayref
    346 !          day_ini = dayref
    347 !          itau_dyn = 0
    348 !          itau_phy = 0
    349 !          time_0 = 0.
    350 !          write(lunout,*)
    351 !     .   'GCM: On reinitialise a la date lue dans gcm.def'
    352 !        endif
    353 !      ELSE
    354 !        raz_date = 0
    355 !      endif
    356 
    357257#ifdef CPP_IOIPSL
    358258  mois = 1
    359259  heure = 0.
    360 ! Ce n'est defini pour l'instant que pour la Terre...
    361   if (planet_type.eq.'earth') then
    362     call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
    363     jH_ref = jD_ref - int(jD_ref)
    364     jD_ref = int(jD_ref)
    365 
    366     call ioconf_startdate(INT(jD_ref), jH_ref)
    367 
    368     write(lunout,*)'DEBUG'
    369     write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
    370     write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
    371     call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
    372     write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
    373     write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
    374 !  else if (planet_type.eq.'titan') then
    375 !    jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
    376   else
    377 ! A voir pour Titan et Venus
    378     jD_ref=0
    379     jH_ref=0
    380     write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
    381     write(lunout,*)jD_ref,jH_ref
    382   endif ! planet_type
    383260#else
    384261  ! Ehouarn: we still need to define JD_ref and JH_ref
    … …  
    410287  !
    411288  !-----------------------------------------------------------------------
    412   !   Initialisation de la dissipation :
    413   !   ----------------------------------
    414 
    415 !ED18 CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
    416 !                  tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
    417 
    418   !-----------------------------------------------------------------------
    419289  !   Initialisation des I/O :
    420290  !   ------------------------
    … …  
    434304  WRITE(lunout,'(a,i7,a,i7)') &
    435305               "run from day ",day_ini,"  to day",day_end
    436 
    437 #ifdef CPP_IOIPSL
    438   ! Ce n'est defini pour l'instant que pour la Terre...
    439   if (planet_type.eq.'earth') then
    440     call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
    441     write (lunout,301)jour, mois, an
    442     call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
    443     write (lunout,302)jour, mois, an
    444   else
    445   ! A voir pour Titan et Venus
    446     write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
    447   endif ! planet_type
    448 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
    449 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
    450 #endif
    451 
    452306
    453307  !-----------------------------------------------------------------------
    … …  
    467321
    468322
    469   if (planet_type=="mars") then
    470     ! For Mars we transmit day_ini
    471     CALL dynredem0("restart.nc", day_ini, phis)
    472   else
    473     CALL dynredem0("restart.nc", day_end, phis)
    474   endif
     323  ! Note that for Mars we transmit day_ini
     324  CALL dynredem0("restart.nc", day_end, phis)
    475325  ecripar = .TRUE.
    476326
    … …  
    503353  istphy=istdyn/iphysiq     
    504354
    505 
    506   !
    507355  !-----------------------------------------------------------------------
    508356  !   Integration temporelle du modele :
    509357  !   ----------------------------------
    510358
    511   !       write(78,*) 'ucov',ucov
    512   !       write(78,*) 'vcov',vcov
    513   !       write(78,*) 'teta',teta
    514   !       write(78,*) 'ps',ps
    515   !       write(78,*) 'q',q
    516 
    517359  CALL leapfrog_nogcm(ucov,vcov,teta,ps,masse,phis,q,time_0)
    518360
Note: See TracChangeset for help on using the changeset viewer.