Changeset 3539


Ignore:
Timestamp:
Dec 9, 2024, 3:39:31 PM (6 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
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified 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

  • TabularUnified 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)

  • TabularUnified 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

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90

    r3455 r3539  
    3131!$OMP THREADPRIVATE(newtonian,force_cpp,cpp_mugaz_mode,testradtimes,rayleigh)
    3232      logical,save :: stelbbody
    33       logical,save :: nearn2cond
    3433      logical,save :: tracer
    3534      logical,save :: mass_redistrib
    36 !$OMP THREADPRIVATE(stelbbody,nearn2cond,tracer,mass_redistrib)
     35!$OMP THREADPRIVATE(stelbbody,tracer,mass_redistrib)
    3736      logical,save :: varactive
    3837      logical,save :: varfixed
    … …  
    8584      logical,save :: convergeps,conservn2,condensn2,no_n2frost
    8685!$OMP THREADPRIVATE(convergeps,conservn2,condensn2,no_n2frost)
     86      logical,save :: conservch4
     87!$OMP THREADPRIVATE(conservch4)
    8788      logical,save :: kmix_proffix
    8889!$OMP THREADPRIVATE(kmix_proffix)

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/ch4surf.F

    r3421 r3539  
    1       SUBROUTINE ch4surf(ngrid,nlayer,nq,ptimestep,tsurf,pdtsurf,
    2      &           pplev,pdpsurf,pq,pdq,pqsurf,pdqsurf,pdqch4,pdqsch4)
     1      SUBROUTINE ch4surf(ngrid,nlayer,nq,ptimestep,capcal,tsurf,
     2     &   pdtsurf,pplev,pdpsurf,pq,pdq,pqsurf,pdqsurf,pdqch4,pdqsch4)
    33 
    44      use callkeys_mod, only: dayfrac, thresh_non2
    … …  
    77      use comsaison_h, only: fract
    88      use planete_mod, only: z0
    9       use tracer_h, only: igcm_ch4_gas,igcm_ch4_ice,igcm_n2,mmol
     9      use tracer_h, only: igcm_ch4_gas,igcm_ch4_ice,igcm_n2,mmol,lw_ch4
    1010      IMPLICIT NONE
    1111
    … …  
    2929
    3030      ! input :
     31      REAL capcal(ngrid)           ! surface heat capacity (J m-2 K-1)
    3132      REAL tsurf(ngrid)
    3233      REAL pplev(ngrid,nlayer+1)
    … …  
    6263      cdrag=(vonk/log(alt/z00))**2
    6364     
    64       u=6.  ! 6
    65       v=3.  ! 3
     65      u=0.3  ! 6
     66      v=0.4  ! 3
    6667      uv=sqrt(u**2+v**2)
    6768       
    … …  
    9798     &             mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.-0.6)/0.3
    9899         gamm(ig)=max(gamm(ig),1.)
    99          qsat(ig)=qsat(ig)*gamm(ig)
     100         !qsat(ig)=qsat(ig)*gamm(ig)
    100101
    101102         rho = zpsrf(ig) / (r *  tsurf(ig) )
    … …  
    125126         !! Security to avoid large changes in temperatures due to
    126127         !latent heat                   
    127          if (pdqsch4(ig)*ptimestep.gt.0.25) then
    128             pdqsch4(ig)=0.25/ptimestep
     128         if (lw_ch4*pdqsch4(ig)*ptimestep/capcal(ig).gt.1.) then
     129            pdqsch4(ig)=1./(lw_ch4*ptimestep)*capcal(ig)
    129130         endif
    130          if (pdqsch4(ig)*ptimestep.lt.-0.25) then
    131             pdqsch4(ig)=-0.25/ptimestep
     131         if (lw_ch4*pdqsch4(ig)*ptimestep/capcal(ig).lt.-1.) then
     132            pdqsch4(ig)=-1./(lw_ch4*ptimestep)*capcal(ig)
    132133         endif
     134         !if (pdqsch4(ig)*ptimestep.gt.0.25) then
     135         !   pdqsch4(ig)=0.25/ptimestep
     136         !endif
     137         !if (pdqsch4(ig)*ptimestep.lt.-0.25) then
     138         !   pdqsch4(ig)=-0.25/ptimestep
     139         !endif
    133140
    134141         !! Atm tendency

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90

    r3477 r3539  
    1313  USE callkeys_mod, only: fast,ch4lag,latlag,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag
    1414  USE vertical_layers_mod, ONLY: ap,bp
    15   use geometry_mod, only: latitude
     15  USE geometry_mod, only: latitude, longitude, cell_area
     16  use planetwide_mod, only: planetwide_sumval
    1617
    1718
    … …  
    8788  REAL tcond_n2
    8889  REAL pcond_n2
    89   REAL glob_average2d         ! function
    9090  REAL zqn2(klon,klev) ! N2 MMR used to compute Tcond/zqn2
    9191  REAL ztcond (klon,klev)
    … …  
    123123
    124124  REAL zplevhist(klon)
     125  REAL zplevhistglob
    125126  REAL zplevnew(klon)
    126127  REAL zplev(klon)
    … …  
    157158  INTEGER,SAVE :: i_n2ice=0       ! n2 ice
    158159  CHARACTER(LEN=20) :: tracername ! to temporarily store text
    159   logical olkin  ! option to prevent N2 ice effect in the south
    160   DATA olkin/.false./
    161   save olkin
    162160
    163161  CHARACTER(len=10) :: tname
    … …  
    176174     ENDDO
    177175     IF (fast) THEN
    178         p00=glob_average2d(kp) ! mean pres at ref level
     176        ! mean pres at ref level
     177        call planetwide_sumval(kp(:)*cell_area(:)/totarea_planet,p00)
    179178     ENDIF
    180179
    … …  
    297296     nsubtimestep= 1 ! if more then kp and p00 have to be calculated
    298297                     ! for nofast mode
    299   ENDIF
     298  ENDIF ! fast
    300299  subtimestep=ptimestep/float(nsubtimestep)
    301300
    … …  
    324323       ENDDO
    325324       ! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code)
    326        globzplevnew=glob_average2d(zplevnew)
     325       call planetwide_sumval(zplevnew(:)*cell_area(:)/totarea_planet,globzplevnew)
    327326       DO ig=1,klon
    328327         zplev(ig)=kp(ig)*globzplevnew/p00  ! 2)
    … …  
    335334          zplevhist(ig)=zplev(ig)
    336335       ENDDO
    337        zplev=zplev*globzplevnew/glob_average2d(zplevhist)   ! 4)
     336       call planetwide_sumval(zplevhist(:)*cell_area(:)/totarea_planet,zplevhistglob)
     337       zplev=zplev*globzplevnew/zplevhistglob   ! 4)
    338338       DO ig=1,klon    ! 5)
    339339         zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))*(ptsrf(ig)**4-ztsrf(ig)**4)
    … …  
    343343     ENDIF  ! (itsub=1)
    344344
    345    DO ig=1,klon
    346      ! forecast of frost temperature ztcondsol
    347      !ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
    348      ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
    349 
    350 !    Loop over where we have condensation / sublimation
    351      IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &    ! ground cond
     345     DO ig=1,klon
     346       ! forecast of frost temperature ztcondsol
     347       !ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
     348       ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
     349
     350       ! Loop over where we have condensation / sublimation
     351       IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &    ! ground cond
    352352                ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &    ! ground sublim
    353353               (zpicen2(ig) .GT. 0.))) THEN
    354         condsub(ig) = .true.    ! condensation or sublimation
    355 
    356 !     Condensation or partial sublimation of N2 ice
    357         if (ztsrf(ig) .LT. ztcondsol(ig)) then  ! condensation
    358 !             Include a correction to account for the cooling of air near
    359 !             the surface before condensing:
    360          zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
    361          /((lw_n2+cpp*(zt(ig,1)-ztcondsol(ig)))*subtimestep)
    362         else                                    ! sublimation
     354         condsub(ig) = .true.    ! condensation or sublimation
     355
     356!        Condensation or partial sublimation of N2 ice
     357         if (ztsrf(ig) .LT. ztcondsol(ig)) then  ! condensation
     358!          Include a correction to account for the cooling of air near
     359!          the surface before condensing:
     360           if (fast) then
     361             zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
     362             /(lw_n2*subtimestep)
     363           else
     364             zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
     365             /((lw_n2+cpp*(zt(ig,1)-ztcondsol(ig)))*subtimestep)
     366           endif
     367         else                                    ! sublimation
    363368          zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
    364369          /(lw_n2*subtimestep)
    365         end if
    366 
    367         pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
    368 
    369 !     partial sublimation of N2 ice
    370 !     If the entire N_2 ice layer sublimes
    371 !     (including what has just condensed in the atmosphere)
    372         IF((zpicen2(ig)/subtimestep).LE. &
     370         end if
     371
     372         pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
     373
     374!        partial sublimation of N2 ice
     375!        If the entire N_2 ice layer sublimes
     376!        (including what has just condensed in the atmosphere)
     377         IF((zpicen2(ig)/subtimestep).LE. &
    373378            -zcondices(ig))THEN
    374379           zcondices(ig) = -zpicen2(ig)/subtimestep
    375380           pdtsrfc(ig)=(lw_n2/pcapcal(ig))*       &
    376381               (zcondices(ig))
    377         END IF
    378 
    379 !     Changing N2 ice amount and pressure
    380 
    381         pdicen2(ig) = zcondices(ig)
    382         pdpsrf(ig)   = -pdicen2(ig)*g
    383     !    pdpsrf(ig)   = 0. ! OPTION to check impact N2 sub/cond
    384         IF (fast.and.(zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001)) then
     382         END IF
     383
     384!        Changing N2 ice amount and pressure
     385         pdicen2(ig) = zcondices(ig)
     386         pdpsrf(ig)   = -pdicen2(ig)*g
     387         ! pdpsrf(ig)   = 0. ! OPTION to check impact N2 sub/cond
     388         IF (fast.and.(zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001)) then
    385389            pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
    386390            pdicen2(ig)=-pdpsrf(ig)/g
    387         ENDIF
    388 
    389      ELSE  ! no condsub
     391         ENDIF
     392
     393       ELSE  ! no condsub
     394
     395        !IF (zpicen2(ig) .GT. 0.) THEN
     396        ! print*,'TB24 no condsub',latitude(ig)*180./pi,longitude(ig)*180./pi
     397        !ENDIF
    390398        pdpsrf(ig)=0.
    391399        pdicen2(ig)=0.
    392400        pdtsrfc(ig)=0.
    393      ENDIF
    394    ENDDO ! ig
    395   enddo ! subtimestep
     401       ENDIF
     402     ENDDO ! ig
     403  enddo ! itsub,nsubtimestep
    396404
    397405!     Updating pressure, temperature and ice reservoir
    … …  
    403411
    404412     pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
     413
     414     !IF (picen2(ig) .GT. 0. .and. abs(pdicen2(ig)).lt.1.e-12) THEN
     415     !    print*,'TB24 low pdq',latitude(ig)*180./pi,longitude(ig)*180./pi
     416     !ENDIF
     417
     418
    405419
    406420! security
    … …  
    413427
    414428     if(.not.picen2(ig).ge.0.) THEN
    415 !           if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
     429!       if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
    416430          print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep
    417 !              pdicen2(ig)= -picen2(ig)/ptimestep
    418 !           else
     431!         pdicen2(ig)= -picen2(ig)/ptimestep
     432!       else
    419433          picen2(ig)=0.0
    420 !           endif
     434!       endif
    421435     endif
    422   ENDDO
     436  ENDDO ! ig
    423437
    424438! ***************************************************************
    … …  
    433447!  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
    434448!  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
    435         zmflux(1) = -zcondices(ig)
    436         DO l=1,klev
     449     zmflux(1) = -zcondices(ig)
     450     DO l=1,klev
    437451         zmflux(l+1) = zmflux(l) -zcondicea(ig,l)  &
    438452         + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
    439453! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
    440454      if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
    441         END DO
     455     END DO
    442456
    443457! Mass of each layer
    444458! ------------------
    445        DO l=1,klev
    446          masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
    447        END DO
    448 
     459     DO l=1,klev
     460        masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
     461     END DO
    449462
    450463!  Corresponding fluxes for T,U,V,Q
    … …  
    456469!           the other tendencies
    457470!           Estimation of subtimestep (using only  the first layer, the most critical)
    458         dtmax=ptimestep
    459         if (zmflux(1).gt.1.e-20) then
    460             dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
    461         endif
    462         nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
    463         subtimestep=ptimestep/float(nsubtimestep)
    464 
    465 !          New flux for each subtimestep
    466        do l=1,klev+1
    467             w(l)=-zmflux(l)*subtimestep
    468        enddo
    469 !          initializing variables that will vary during subtimestep:
    470        do l=1,klev
     471     dtmax=ptimestep
     472     if (zmflux(1).gt.1.e-20) then
     473         dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
     474     endif
     475     nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
     476     subtimestep=ptimestep/float(nsubtimestep)
     477
     478!    New flux for each subtimestep
     479     do l=1,klev+1
     480         w(l)=-zmflux(l)*subtimestep
     481     enddo
     482!    initializing variables that will vary during subtimestep:
     483     do l=1,klev
    471484           ztc(l)  =pt(ig,l)
    472485           zu(l)   =pu(ig,l)
    … …  
    475488              zq(l,iq) = pq(ig,l,iq)
    476489           enddo
    477        end do
    478 
    479 !          loop over nsubtimestep
    480 !          """"""""""""""""""""""
    481        do itsub=1,nsubtimestep
    482 !             Progressively adding tendancies from other processes.
    483           do l=1,klev
     490     end do
     491
     492!    loop over nsubtimestep
     493!    """"""""""""""""""""""
     494     do itsub=1,nsubtimestep
     495!        Progressively adding tendancies from other processes.
     496         do l=1,klev
    484497             ztc(l)  =ztc(l)   +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
    485498             zu(l)   =zu(l)   +pdu( ig,l) * subtimestep
    … …  
    488501                zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
    489502             enddo
    490            end do
    491 
    492 !             Change of mass in each layer
    493           do l=1,klev
     503         end do
     504
     505!        Change of mass in each layer
     506         do l=1,klev
    494507             masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))&
    495508                                             /(g*pplev(ig,1))
    496           end do
    497 
    498 !             Value transfert at the surface interface when condensation sublimation:
    499 
    500           if (zmflux(1).lt.0) then
    501 !               Surface condensation
    502             zum(1)= zu(1)
    503             zvm(1)= zv(1)
    504             ztm(1) = ztc(1)
    505           else
    506 !               Surface sublimation:
    507             ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
    508             zum(1) = 0
    509             zvm(1) = 0
    510           end if
    511           do iq=1,nq
    512               zqm(1,iq)=0. ! most tracer do not condense !
    513           enddo
    514 !             Special case if the tracer is n2 gas
    515           if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
    516 
    517            ztm(2:klev+1)=0.
    518            zum(2:klev+1)=0.
    519            zvm(2:klev+1)=0.
    520            zqm1(1:klev+1)=0.
    521 
    522 !             Van Leer scheme:
    523           call vl1d(klev,ztc,2.,masse,w,ztm)
    524           call vl1d(klev,zu ,2.,masse,w,zum)
    525           call vl1d(klev,zv ,2.,masse,w,zvm)
     509         end do
     510
     511         ztm(2:klev+1)=0.
     512         zum(2:klev+1)=0.
     513         zvm(2:klev+1)=0.
     514         zqm1(1:klev+1)=0.
     515
     516!        Van Leer scheme:
     517         call vl1d(klev,ztc,2.,masse,w,ztm)
     518         call vl1d(klev,zu ,2.,masse,w,zum)
     519         call vl1d(klev,zv ,2.,masse,w,zvm)
    526520         do iq=1,nq
    527521           do l=1,klev
    … …  
    533527            zqm(l,iq)=zqm1(l)
    534528           enddo
    535           enddo
    536 
    537 !             Correction to prevent negative value for qn2
    538           if (igcm_n2.ne.0) then
     529         enddo
     530
     531!        Correction to prevent negative value for qn2
     532         if (igcm_n2.ne.0) then
    539533            zqm(1,igcm_n2)=1.
    540534            do l=1,klev-1
    … …  
    566560          if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
    567561
    568            !!! Source haze: 0.02 pourcent when n2 sublimes
     562          !!! Source haze: 0.02 pourcent when n2 sublimes
    569563          IF (source_haze) THEN
    570564               IF (pdicen2(ig).lt.0) THEN
    … …  
    611605          END DO
    612606
    613 !             Tendencies on Q
     607!         Tendencies on Q
    614608          do iq=1,nq
    615609            if (iq.eq.igcm_n2) then
    616 !                 SPECIAL Case when the tracer IS N2 :
     610!             SPECIAL Case when the tracer IS N2 :
    617611              DO l=1,klev
    618612              pdqc(ig,l,iq)= (1/masse(l)) *        &
    … …  
    630624            end if
    631625          enddo
    632 !             Update variables at the end of each subtimestep.
     626!         Update variables at the end of each subtimestep.
    633627          do l=1,klev
    634628            ztc(l)  =ztc(l)  + zdtsig(ig,l)  *subtimestep
    … …  
    639633            enddo
    640634          end do
    641        enddo   ! loop on nsubtimestep
    642 !          Recomputing Total tendencies
    643        do l=1,klev
     635
     636     enddo   ! loop on nsubtimestep
     637
     638!    Recomputing Total tendencies
     639     do l=1,klev
    644640         pdtc(ig,l)   = (ztc(l) - zt(ig,l) )/ptimestep
    645641         pduc(ig,l)   = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
    … …  
    659655
    660656         enddo
    661        end do
     657     end do
    662658! *******************************TEMPORAIRE ******************
    663        if (klon.eq.1) then
     659     if (klon.eq.1) then
    664660              write(*,*) 'nsubtimestep=' ,nsubtimestep
    665661           write(*,*) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
    … …  
    671667              write(*,*) 'dq*Dt l=1' ,  pdq(1,1,1)*ptimestep
    672668              write(*,*) 'dqcond*Dt l=1' ,  pdqc(1,1,1)*ptimestep
    673        end if
     669     end if
    674670
    675671! ***********************************************************
    … …  
    677673   END DO  ! loop on ig
    678674  endif ! not fast
    679 
    680 ! ************ Option Olkin to prevent N2 effect in the south********
    681 112   continue
    682   if (olkin) then
    683   DO ig=1,klon
    684      if (latitude(ig).lt.0) then
    685        pdtsrfc(ig) = max(0.,pdtsrfc(ig))
    686        pdpsrf(ig) = 0.
    687        pdicen2(ig)  = 0.
    688        do l=1,klev
    689          pdtc(ig,l)  =  max(0.,zdtlatent(ig,l))
    690          pduc(ig,l)  = 0.
    691          pdvc(ig,l)  = 0.
    692          do iq=1,nq
    693            pdqc(ig,l,iq) = 0
    694          enddo
    695        end do
    696      end if
    697   END DO
    698   end if
    699 ! *******************************************************************
    700675
    701676! ***************************************************************
    … …  
    764739!-------------------------------------------------------------------------
    765740
    766   real function glob_average2d(var)
    767 !   Calculates the global average of variable var
    768   use comgeomfi_h
    769   use dimphy, only: klon
    770   USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
    771   use geometry_mod, only: cell_area, latitude
    772 
    773   implicit none
    774 
    775 !      INTEGER klon
    776   REAL var(klon)
    777   INTEGER ig
    778 
    779   glob_average2d = 0.
    780   DO ig=2,klon-1
    781        glob_average2d = glob_average2d + var(ig)*cell_area(ig)
    782   END DO
    783   glob_average2d = glob_average2d + var(1)*cell_area(1)*nbp_lon
    784   glob_average2d = glob_average2d + var(klon)*cell_area(klon)*nbp_lon
    785   glob_average2d = glob_average2d/(totarea+(cell_area(1)+cell_area(klon))*(nbp_lon-1))
    786 
    787   end function glob_average2d
    788 
    789 ! *****************************************************************
    790 
    791741  subroutine vl1d(klev,q,pente_max,masse,w,qm)
    792742!
    … …  
    835785  enddo
    836786
    837      dzq(1)=0.
    838      dzq(klev)=0.
    839 
    840    do l = 1,klev-1
     787  dzq(1)=0.
     788  dzq(klev)=0.
     789
     790  do l = 1,klev-1
    841791
    842792!         Regular scheme (transfered mass < layer mass)
    … …  
    895845         end if
    896846      end if
    897    enddo
     847  enddo
    898848
    899849   return

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/cosurf.F

    r3421 r3539  
    5959      cdrag=(vonk/log(alt/z00))**2
    6060     
    61       u=6.
    62       v=3.
     61      u=0.3
     62      v=0.4
    6363      uv=sqrt(u**2+v**2)
    6464       

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90

    r3477 r3539  
    2121  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
    2222                          mugaz, pi, avocado
    23   use planete_mod, only: nres
    2423  use planetwide_mod, only: planetwide_sumval
    2524  use callkeys_mod
    … …  
    203202     call getin_p("convergeps",convergeps)
    204203     if (is_master) write(*,*) trim(rname)//" convergeps = ",convergeps
     204
    205205     conservn2=.false. ! default value
    206206     call getin_p("conservn2",conservn2)
    207207     if (is_master) write(*,*) trim(rname)//" conservn2 = ",conservn2
     208     conservch4=.false. ! default value
     209     call getin_p("conservch4",conservch4)
     210     if (is_master) write(*,*) trim(rname)//" conservch4 = ",conservch4
    208211
    209212     if (is_master) write(*,*) trim(rname)//"KBO runs (eris, makemake) ?"
    … …  
    227230     call getin_p("noseason_day",noseason_day)
    228231     if (is_master) write(*,*) trim(rname)//": noseason_day=", noseason_day
    229 
    230      if (is_master) write(*,*) trim(rname)//": Tidal resonance ratio ?"
    231      nres=0          ! default value
    232      call getin_p("nres",nres)
    233      if (is_master) write(*,*) trim(rname)//": nres = ",nres
    234232
    235233     if (is_master) write(*,*) trim(rname)//&
    … …  
    310308     if (is_master) write(*,*) trim(rname)//&
    311309       ": call turbulent vertical diffusion ?"
    312      calldifv=.true. ! default value
     310     calldifv=.false. ! default value
    313311     call getin_p("calldifv",calldifv)
    314312     if (is_master) write(*,*) trim(rname)//": calldifv = ",calldifv
    315313
    316314     if (is_master) write(*,*) trim(rname)//": use turbdiff instead of vdifc ?"
    317      UseTurbDiff=.true. ! default value
     315     UseTurbDiff=.false. ! default value
    318316     call getin_p("UseTurbDiff",UseTurbDiff)
    319317     if (is_master) write(*,*) trim(rname)//": UseTurbDiff = ",UseTurbDiff
    320318
    321319     if (is_master) write(*,*) trim(rname)//": call convective adjustment ?"
    322      calladj=.true. ! default value
     320     calladj=.false. ! default value
    323321     call getin_p("calladj",calladj)
    324322     if (is_master) write(*,*) trim(rname)//": calladj = ",calladj
    … …  
    14771475         call abort_physic(rname, 'for now, haze/rad proffix only works w aerohaze=T', 1)
    14781476     endif
    1479       if (condcosurf.and.no_n2frost) then
     1477      if (carbox.and.condcosurf.and.no_n2frost) then
    14801478        call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
    14811479      end if

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/iniwritesoil.F90

    r3184 r3539  
    195195text="Time"
    196196ierr=NF_PUT_ATT_TEXT(nid,varid,"long_name",len_trim(text),text)
    197 text="days since 0000-01-01 00:00:00"
     197text="days since 0000-00-00 00:00:00"
    198198ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
    199199

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/phyredem.F90

    r3184 r3539  
    77subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
    88                         phystep,day_ini,time,airefi, &
    9                          alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
     9                         alb,pzmea,pzstd,pzsig,pzgam,pzthe)
    1010! create physics restart file and write time-independent variables
    1111  use comsoil_h, only: volcapa, mlayer
    … …  
    3636  real,intent(in) :: airefi(ngrid)
    3737  real,intent(in) :: alb(ngrid)
    38   real,intent(in) :: ith(ngrid,nsoil)
    3938  real,intent(in) :: pzmea(ngrid)
    4039  real,intent(in) :: pzstd(ngrid)
    … …  
    125124  call put_field("ZTHE","Relief: theta parameter",zthe)
    126125 
    127   ! Soil thermal inertia
    128   call put_field("inertiedat","Soil thermal inertia",ith)
    129  
    130126  ! Close file
    131127  call close_restartphy
    … …  
    134130
    135131subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
    136                     phystep,time,tsurf,tsoil,emis,q2,qsurf)
     132                    phystep,time,tsurf,tsoil,inertiesoil, &
     133                    emis,q2,qsurf)
    137134  ! write time-dependent variable to restart file
    138135  use iostart, only : open_restartphy, close_restartphy, &
    … …  
    152149  real,intent(in) :: tsoil(ngrid,nsoil)
    153150  real,intent(in) :: emis(ngrid)
     151  real,intent(in) :: inertiesoil(ngrid,nsoil)
    154152  real,intent(in) :: q2(ngrid,nlay+1)
    155153  real,intent(in) :: qsurf(ngrid,nq)
    … …  
    166164  ! Surface temperature
    167165  call put_field("tsurf","Surface temperature",tsurf)
     166
     167  ! Soil inertia
     168  call put_field("inertiedat","Soil thermal inertia",inertiesoil)
    168169 
    169170  ! Soil temperature

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/phys_state_var_mod.F90

    r3195 r3539  
    8181
    8282
    83       real,allocatable,dimension(:,:),save :: qsurf_hist
    8483      real,allocatable,dimension(:,:,:),save :: nueffrad ! Aerosol effective radius variance. By RW
    8584      real,allocatable,dimension(:,:,:),save :: reffrad
    86 !$OMP THREADPRIVATE(qsurf_hist,nueffrad,reffrad)
    8785
    8886      real,dimension(:,:),allocatable,save :: dEzdiff   ! Turbulent diffusion heating (W.m-2)
    … …  
    143141        ALLOCATE(cloudfrac(klon,klev))
    144142        ALLOCATE(totcloudfrac(klon))
    145         ALLOCATE(qsurf_hist(klon,nqtot))
    146143        ALLOCATE(reffrad(klon,klev,naerkind))
    147144        ALLOCATE(nueffrad(klon,klev,naerkind))
    … …  
    221218        DEALLOCATE(cloudfrac)
    222219        DEALLOCATE(totcloudfrac)
    223         DEALLOCATE(qsurf_hist)
    224220        DEALLOCATE(reffrad)
    225221        DEALLOCATE(nueffrad)

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

    r3522 r3539  
    2929      use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
    3030      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
    31       use geometry_mod, only: latitude, longitude, cell_area
     31      use geometry_mod, only: latitude, longitude, cell_area, &
     32                          cell_area_for_lonlat_outputs
    3233      USE comgeomfi_h, only: totarea, totarea_planet
    3334      USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
    … …  
    4344      use mod_phys_lmdz_para, only : is_master
    4445      use planete_mod, only: apoastr, periastr, year_day, peri_day, &
    45                             obliquit, nres, z0, adjust, tpal
     46                            obliquit, z0, adjust, tpal
    4647      use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
    4748      use time_phylmdz_mod, only: daysec
    … …  
    5455                              lwrite, mass_redistrib, meanOLR,                &
    5556                              fast,fasthaze,haze,metcloud,monoxcloud,         &
    56                               n2cond,nearn2cond,noseason_day,conservn2,       &
     57                              n2cond,noseason_day,conservn2,conservch4,       &
    5758                              convergeps,kbo,triton,paleo,paleoyears,glaflow, &
    5859                              carbox, methane,condmetsurf,condcosurf,         &
    … …  
    7879      use condensation_generic_mod, only: condensation_generic
    7980      use datafile_mod, only: datadir
    80 #ifndef MESOSCALE
    8181      USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
    8282      use mod_phys_lmdz_omp_data, ONLY: is_omp_master
    83 #else
    84       use comm_wrf, only : comm_HR_SW, comm_HR_LW,  &
    85          comm_CLOUDFRAC,comm_TOTCLOUDFRAC, comm_RH, &
    86          comm_HR_DYN,     &
    87          comm_DQICE,comm_DQVAP,comm_ALBEQ,          &
    88          comm_FLUXTOP_DN,comm_FLUXABS_SW,           &
    89          comm_FLUXTOP_LW,comm_FLUXSURF_SW,          &
    90          comm_FLUXSURF_LW,comm_FLXGRD,              &
    91          comm_DTRAIN,comm_DTLSC,                    &
    92          comm_LATENT_HF
    93 
    94 #endif
     83      USE mod_grid_phy_lmdz, ONLY: regular_lonlat, grid_type, unstructured
    9584
    9685#ifdef CPP_XIOS
    … …  
    11099!     -------
    111100!     Central subroutine for all the physics parameterisations in the
    112 !     universal model. Originally adapted from the Mars LMDZ model.
    113 !
    114 !     The model can be run without or with tracer transport
     101!     Pluto model. Originally adapted from the Mars LMDZ model.
     102!
     103!     The model can be run with 1 (N2) or more tracer transport
    115104!     depending on the value of "tracer" in file "callphys.def".
    116 !
    117105!
    118106!   It includes:
    … …  
    122110!         I.2 Initialization for every call to physiq.
    123111!
    124 !      II. Compute radiative transfer tendencies (longwave and shortwave) :
     112!      II.1 Thermosphere
     113!      II.2 Compute radiative transfer tendencies (longwave and shortwave) :
    125114!         II.a Option 1 : Call correlated-k radiative transfer scheme.
    126 !         II.b Option 2 : Call Newtonian cooling scheme.
    127 !         II.! Option 3 : Atmosphere has no radiative effect.
    128 !
    129 !      III. Vertical diffusion (turbulent mixing) :
     115!         II.b Option 2 : Atmosphere has no radiative effect.
     116!
     117!      III. Vertical diffusion (turbulent mixing)
    130118!
    131119!      IV. Convection :
    132 !         IV.a Thermal plume model
    133 !         IV.b Dry convective adjusment
    134 !
    135 !      V. Condensation and sublimation of gases (currently just N2).
     120!         IV.a Dry convective adjusment
     121!
     122!      V. Condensation and sublimation of gases.
    136123!
    137124!      VI. Tracers
    138 !         VI.1. Water and water ice.
    139 !         VI.2  Photochemistry
    140 !         VI.3. Aerosols and particles.
    141 !         VI.4. Updates (pressure variations, surface budget).
    142 !         VI.5. Slab Ocean.
    143 !         VI.6. Surface Tracer Update.
     125!         VI.1. Aerosols and particles.
     126!         VI.2. Updates (pressure variations, surface budget).
     127!         VI.3. Surface Tracer Update.
    144128!
    145129!      VII. Surface and sub-surface soil temperature.
    … …  
    210194!           Photochemical core developped by F. Lefevre: B. Charnay (2017)
    211195!           Purge for Pluto model : A. Falco (2024)
     196!           Adapting to Pluto : A. Falco, T. Bertrand (2024)
    212197!==================================================================
    213198
    … …  
    264249      REAL,save  :: tcond1p4Pa
    265250      DATA tcond1p4Pa/38/
    266 
    267 
    268251
    269252! Local variables :
    … …  
    297280      REAL,SAVE :: glastep=20   ! step in pluto day to spread glacier
    298281
    299 
    300282!     Aerosol (dust or ice) extinction optical depth  at reference wavelength
    301283!     for the "naerkind" optically active aerosols:
    … …  
    318300      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
    319301
    320 
    321 ! VARIABLES for the thermal plume model (AF24: deleted)
    322 
    323302! TENDENCIES due to various processes :
    324303
    … …  
    328307      real zdtsurfc(ngrid)    ! Condense_n2 routine.
    329308      real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
    330       ! real zdtsurf_hyd(ngrid) ! Hydrol routine.
    331309
    332310      ! For Atmospheric Temperatures : (K/s)
    … …  
    359337      REAL,allocatable,save :: zdqschim(:,:)  ! Calchim_asis routine
    360338!$OMP THREADPRIVATE(zdqchim,zdqschim)
    361 
    362339
    363340      !! PLUTO variables
    … …  
    400377      REAL zfluxuv                     ! Lyman alpha flux at 1AU
    401378
    402 
    403 
    404379      REAL array_zero1(ngrid)
    405380      REAL array_zero2(ngrid,nlayer)
    … …  
    417392      real zdmassmr_col(ngrid)    ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
    418393      real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
    419 
    420 
    421394
    422395! Local variables for LOCAL CALCULATIONS:
    … …  
    431404      real gmplanet
    432405      real taux(ngrid),tauy(ngrid)
    433 
    434406
    435407! local variables for DIAGNOSTICS : (diagfi & stat)
    … …  
    471443      real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
    472444
    473 
    474445!     to test energy conservation (RW+JL)
    475446      real mass(ngrid,nlayer),massarea(ngrid,nlayer)
    … …  
    479450
    480451!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
    481 
    482452      real dtmoist_max,dtmoist_min
    483453      real dItot, dItot_tmp, dVtot, dVtot_tmp
    484 
    485 
    486454      real dWtot, dWtot_tmp, dWtots, dWtots_tmp
    487 
    488455      real psat_tmp ! AF24: to remove?
    489456
    … …  
    512479      REAL zustrhi(ngrid), zvstrhi(ngrid)
    513480
    514 
    515481      real :: tsurf2(ngrid)
    516482!!      real :: flux_o(ngrid),flux_g(ngrid)
    … …  
    518484      real :: flux_sens_lat(ngrid)
    519485      real :: qsurfint(ngrid,nq)
    520 #ifdef MESOSCALE
    521       REAL :: lsf_dt(nlayer)
    522       REAL :: lsf_dq(nlayer)
    523 #endif
    524486
    525487      ! local variables for skin depth check
     488      real :: therm_inertia(ngrid,nsoilmx)
    526489      real :: inertia_min,inertia_max
    527490      real :: diurnal_skin ! diurnal skin depth (m)
    … …  
    556519         ! Allocate saved arrays (except for 1D model, where this has already
    557520         ! been done)
    558 #ifndef MESOSCALE
    559521         if (ngrid>1) call phys_state_var_init(nq)
    560 #endif
    561522
    562523!        Variables set to 0
    … …  
    591552!        Read 'startfi.nc' file.
    592553!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    593 #ifndef MESOSCALE
    594554         call phyetat0(startphy_file,                                 &
    595555                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
    596556                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat)
    597557
    598 #else
    599 
    600          day_ini = pday
    601 #endif
    602 
    603 #ifndef MESOSCALE
    604558         if (.not.startphy_file) then
    605559           ! starting without startfi.nc and with callsoil
    … …  
    611565              alpha=2
    612566              do iloop=0,nsoilmx-1
    613               mlayer(iloop)=lay1*(alpha**(iloop-0.5))
    614           enddo
     567                 mlayer(iloop)=lay1*(alpha**(iloop-0.5))
     568              enddo
    615569              lay1=sqrt(mlayer(0)*mlayer(1))
    616570              alpha=mlayer(1)/mlayer(0)
    … …  
    629583           day_ini=pday
    630584         endif
    631 #endif
     585
    632586         if (pday.ne.day_ini) then
    633587             write(*,*) "ERROR in physiq.F90:"
    … …  
    674628         icount=1
    675629
    676 !        Initialize surface history variable.
    677 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    678          qsurf_hist(:,:)=qsurf(:,:)
    679 
    680 !!         call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
    681 
    682630!        Initialize variable for dynamical heating and zonal wind tendency diagnostic
    683631!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    685633         zuprevious(:,:)=pu(:,:)
    686634
    687 !        Set temperature just above condensation temperature (for Early Mars)
    688 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    689          if(nearn2cond) then
    690             write(*,*)' WARNING! Starting at Tcond+1K'
    691             do l=1, nlayer
    692                do ig=1,ngrid
    693                   pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4     &
    694                       -pt(ig,l)) / ptimestep
    695                enddo
    696             enddo
    697          endif
    698 
    699635         if(meanOLR)then
    700636            call system('rm -f rad_bal.out') ! to record global radiative balance.
    … …  
    703639         endif
    704640
    705 
    706          !! Initialize variables for water cycle ! AF24: removed
    707 
    708641!        Set metallicity for GCS
    709642!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    714647!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    715648
    716 #ifndef MESOSCALE
    717649         if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
    718650            call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
    719651                         ptimestep,pday+nday,time_phys,cell_area,          &
    720                          albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
    721          endif
    722 
    723 !!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1)
    724 #endif
     652                         albedo_bareground,zmea,zstd,zsig,zgam,zthe)
     653         endif
     654
    725655         if (corrk) then
    726656            ! We initialise the spectral grid here instead of
    … …  
    779709      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    780710
    781       if ( .not.nearn2cond ) then
    782          pdt(1:ngrid,1:nlayer) = 0.0
    783       endif
     711      pdt(1:ngrid,1:nlayer) = 0.0
    784712      zdtsurf(1:ngrid)      = 0.0
    785713      pdq(1:ngrid,1:nlayer,1:nq) = 0.0
    … …  
    808736      end if
    809737
    810 
    811738!     Get Lyman alpha flux at specific Ls
    812739      if (haze) then
    … …  
    815742      end if
    816743
    817 
    818744      IF (triton) then
    819745         CALL orbitetriton(zls,zday,dist_star,declin)
    … …  
    822748      ENDIF
    823749
    824 
    825750      if (diurnal) then
    826751              zlss=-2.*pi*(zday-.5)
    … …  
    828753              zlss=9999.
    829754      endif
    830 
    831755
    832756      glat(:) = g !AF24: removed oblateness
    … …  
    892816      endif
    893817
    894 
    895 
    896818! --------------------------------------------------------
    897 !    1.3 thermosphere
     819!    II.1 Thermosphere
    898820! --------------------------------------------------------
    899821
    … …  
    944866
    945867      if (conservn2) then
    946          write(*,*) 'conservn2 thermo'
     868         write(*,*) 'conservn2 thermosphere'
    947869         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
    948870      endif
    … …  
    950872
    951873!---------------------------------
    952 ! II. Compute radiative tendencies
     874! II.2 Compute radiative tendencies
    953875!---------------------------------
    954876!     Saving qsurf to compute paleo flux condensation/sublimation
    … …  
    999921            .and.zday.gt.saveday+1000.   &
    1000922            .and.declin.lt.savedeclin) then
    1001                call globalaverage2d(ngrid,pplev(:,1),globave)
     923               call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
    1002924               if (globave.gt.1.5) then
    1003925                     adjust=adjust+0.005
    … …  
    1012934            .and.zday.gt.saveday+10000.  &
    1013935            .and.declin.gt.savedeclin) then
    1014                call globalaverage2d(ngrid,pplev(:,1),globave)
     936               call planetwide_sumval(pplev(:,1)*cell_area(:)/totarea_planet,globave)
    1015937               if (globave.gt.1.2) then
    1016938                     adjust=adjust+0.005
    … …  
    1026948      call surfprop(ngrid,nq,fract,qsurf,tsurf,        &
    1027949      capcal,adjust,dist_star,flusurfold,ptimestep,zls,&
    1028       albedo,emis,inertiedat)
     950      albedo,emis,therm_inertia)
    1029951      ! do i=2,ngrid
    1030952      !    albedo(i,:) = albedo(1,:)
    … …  
    1032954
    1033955      if (firstcall.and.callsoil) then
    1034          ! AF24 Originally in soil.F, but inertiedat is modified by surfprop
     956         ! AF24 Originally in soil.F, but therm_inertia is modified by surfprop
    1035957         ! Additional checks: is the vertical discretization sufficient
    1036958         ! to resolve diurnal and annual waves?
    1037959         do ig=1,ngrid
    1038960            ! extreme inertia for this column
    1039             inertia_min=minval(inertiedat(ig,:))
    1040             inertia_max=maxval(inertiedat(ig,:))
     961            inertia_min=minval(therm_inertia(ig,:))
     962            inertia_max=maxval(therm_inertia(ig,:))
    1041963            ! diurnal and annual skin depth
    1042964            diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
    … …  
    11701092               endif
    11711093
    1172 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1173 ! II.b Call Newtonian cooling scheme !AF24: removed
    1174 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11751094            else
    11761095! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1177 ! II.! Atmosphere has no radiative effect
     1096! II.b Atmosphere has no radiative effect
    11781097! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11791098               fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
    … …  
    12641183            acond=r/lw_n2
    12651184
    1266             !------------------------------------------------------------------
    1267             ! test methane conservation
    1268             !         if(methane)then
    1269             !           call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
    1270             !     &          ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
    1271             !         endif  ! methane
    1272             !------------------------------------------------------------------
    1273             ! test CO conservation
    1274             !         if(carbox)then
    1275             !           call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
    1276             !     &          ptimestep,pplev,zdqdif,zdqsdif,'CO ',' vdifc ')
    1277             !         endif  ! carbox
    1278             !------------------------------------------------------------------
    1279 
    12801185         ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
    12811186         else if (UseTurbDiff) then
    … …  
    13611266         endif ! end of 'enertest'
    13621267
    1363 
    1364          ! ! Test water conservation. !AF24: removed
    1365 
    13661268      else ! calldifv
    13671269
     1270         ztim1=4.*sigma*ptimestep
     1271         DO ig=1,ngrid
     1272           ztim2=ztim1*emis(ig)*tsurf(ig)**3
     1273           z1=capcal(ig)*tsurf(ig)+  &
     1274           ztim2*tsurf(ig)+ (fluxrad(ig)+fluxgrd(ig))*ptimestep
     1275           z2= capcal(ig)+ztim2
     1276           zdtsurf(ig)=(z1/z2 - tsurf(ig))/ptimestep
     1277
     1278           ! for output:
     1279           !dplanck(ig)=4.*stephan*ptimestep*emis(ig)*tsurf(ig)**3
     1280         ENDDO
     1281
    13681282         ! if(.not.newtonian)then
    1369          zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
     1283         !zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
    13701284
    13711285!        ------------------------------------------------------------------
    … …  
    13741288         if ((methane).and.(fast).and.condmetsurf) THEN
    13751289
    1376             call ch4surf(ngrid,nlayer,nq,ptimestep, &
     1290            call ch4surf(ngrid,nlayer,nq,ptimestep,capcal, &
    13771291               tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, &
    13781292               zdqch4fast,zdqsch4fast)
    … …  
    14041318
    14051319      if (conservn2) then
    1406         write(*,*) 'conservn2 diff'
     1320        write(*,*) 'conservn2 calldifv'
    14071321        call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ &
    14081322                                       dqsurf(:,1)*ptimestep)
     1323      endif
     1324      if (methane.and.conservch4) then
     1325        write(*,*) 'conservch4 calldifv'
     1326        if (fast) then
     1327             call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
     1328                   qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
     1329                   ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' vdifc ')
     1330        else
     1331             call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
     1332                   igcm_ch4_gas,igcm_ch4_ice, &
     1333                   ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
     1334        endif
    14091335      endif
    14101336
    … …  
    14131339!-------------------
    14141340
    1415 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
    1416 ! IV.a Thermal plume model : AF24: removed
    1417 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
    1418 
    14191341! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1420 ! IV.b Dry convective adjustment :
     1342! IV.a Dry convective adjustment :
    14211343! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    14221344
    … …  
    14541376
    14551377      endif ! end of 'calladj'
    1456 !----------------------------------------------
    1457 !   Non orographic Gravity Waves: AF24: removed
    1458 !---------------------------------------------
    14591378
    14601379!-----------------------------------------------
    … …  
    15061425           pdpsrf(:)*ptimestep,qsurf(:,1)+dqsurf(:,1)*ptimestep)
    15071426      endif
     1427      if (methane.and.conservch4) then
     1428        write(*,*) 'conservch4 n2cond'
     1429        if (fast) then
     1430             call testconservfast(ngrid,nlayer,nq,pq(:,1,igcm_ch4_gas),pdq(:,1,igcm_ch4_gas), &
     1431                   qsurf(:,igcm_ch4_ice),dqsurf(:,igcm_ch4_ice), &
     1432                   ptimestep,pplev,zdqch4fast,zdqsch4fast,'CH4',' n2cond')
     1433        else
     1434             call testconserv(ngrid,nlayer,nq,pq,pdq,qsurf,dqsurf, &
     1435                   igcm_ch4_gas,igcm_ch4_ice, &
     1436                   ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond')
     1437        endif
     1438      endif
    15081439
    15091440!---------------------------------------------
    … …  
    15121443
    15131444      if (tracer) then
    1514 
    1515 
    15161445
    15171446!   7a. Methane, CO, and ice
    … …  
    15841513         rice_co(:,:)=0 ! initialization needed for callsedim
    15851514         END IF  ! of IF (carbox)
    1586 
    1587 !------------------------------------------------------------------
    1588 ! test methane conservation
    1589 !         if(methane)then
    1590 !           call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
    1591 !     &         ptimestep,pplev,zdqch4cloud,zdqsch4cloud,'CH4','ch4clou')
    1592 !         endif  ! methane
    1593 !------------------------------------------------------------------
    1594 ! test CO conservation
    1595 !         if(carbox)then
    1596 !           call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
    1597 !     &          ptimestep,pplev,zdqcocloud,zdqscocloud,'CO ','cocloud')
    1598 !         endif  ! carbox
    1599 !------------------------------------------------------------------
    16001515
    16011516!   7b. Haze particle production
    … …  
    16461561      ENDIF
    16471562
    1648 
    1649 
    16501563  ! -------------------------
    16511564  !   VI.3. Aerosol particles
    … …  
    17481661!         call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
    17491662
    1750 !!         call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas))
    1751 !!         call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
    1752 
    1753 
    1754 ! ------------------
    1755 !   VI.5. Slab Ocean !AF24: removed
    1756 ! ------------------
    1757 
    17581663  ! -----------------------------
    17591664  !   VI.6. Surface Tracer Update
    17601665  ! -----------------------------
    17611666
    1762          ! AF24: deleted hydrology
    1763 
    17641667         qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
    1765 
    1766          ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
    1767          ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
    1768          qsurf_hist(:,:) = qsurf(:,:)
    1769 
    1770          ! if(ice_update)then
    1771          !    ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice))
    1772          ! endif
    17731668
    17741669      endif! end of if 'tracer'
    … …  
    18201715!   For diagnostic
    18211716!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1822       DO ig=1,ngrid
     1717      if (.not.fast) then
     1718       DO ig=1,ngrid
    18231719         rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
    18241720         sensiblehf1(ig)=rho(ig,1)*cpp*(0.4/log(zzlay(ig,1)/z0))**2* &
    … …  
    18281724            sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)
    18291725         end if
    1830 
    1831       ENDDO
     1726       ENDDO
     1727      endif
    18321728
    18331729
    … …  
    18501746!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    18511747      if (callsoil) then
    1852          call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
     1748         call soil(ngrid,nsoilmx,.false.,lastcall,therm_inertia,   &
    18531749                   ptimestep,tsurf,tsoil,capcal,fluxgrd)
    18541750      endif
    … …  
    19391835      endif
    19401836
    1941       ! Dynamical heating diagnostic.
    1942       do ig=1,ngrid
    1943          fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp
    1944       enddo
     1837      ! Dynamical heating diagnostic
     1838      fluxdyn(:)=0.0
     1839      if (.not.fast) then
     1840        do ig=1,ngrid
     1841           fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp
     1842        enddo
     1843      endif
    19451844
    19461845      ! Tracers.
    … …  
    19491848      ! Surface pressure.
    19501849      ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
    1951       call globalaverage2d(ngrid,ps,globave)
     1850      call planetwide_sumval(ps(:)*cell_area(:)/totarea_planet,globave)
    19521851
    19531852      ! pressure density !pluto specific
    … …  
    22162115            ! Finally ensure conservation of qsurf
    22172116            DO iq=1,nq
    2218                call globalaverage2d(ngrid,qsurf(:,iq),globaveice(iq))
    2219                call globalaverage2d(ngrid,qsurfpal(:,iq), &
    2220                                              globavenewice(iq))
     2117               call planetwide_sumval(qsurf(:,iq)*cell_area(:)/totarea_planet,globaveice(iq))
     2118               call planetwide_sumval(qsurfpal(:,iq)*cell_area(:)/totarea_planet,globavenewice(iq))
    22212119               IF (globavenewice(iq).gt.0.) THEN
    22222120                  qsurfpal(:,iq)=qsurfpal(:,iq)* &
    … …  
    22632161!              thus we store for time=time+dtvr
    22642162
    2265 #ifndef MESOSCALE
    2266              ! create restartfi
     2163            ! create restartfi
    22672164            if (ngrid.ne.1) then
    22682165               print*, "physdem1pal not yet implemented"
    … …  
    22722169               !      ptimestep,pdaypal, &
    22732170               !      ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, &
    2274                !      cell_area,albedodat,inertiedat,zmea,zstd,zsig, &
     2171               !      cell_area,albedodat,therm_inertia,zmea,zstd,zsig, &
    22752172               !      zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal,  &
    22762173               !      peri_daypal)
    … …  
    22832180               call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
    22842181                           ptimestep,ztime_fin,                    &
    2285                            tsurf,tsoil,emis,q2,qsurf_hist)
     2182                           tsurf,tsoil,therm_inertia,emis,q2,qsurf)
    22862183            endif
    2287 #endif
     2184
    22882185         endif ! end of 'paleo'
    22892186      endif ! end of 'lastcall'
    2290 
    2291 
    2292 ! -----------------------------------------------------------------
    2293 ! WSTATS: Saving statistics
    2294 ! -----------------------------------------------------------------
    2295 ! ("stats" stores and accumulates key variables in file "stats.nc"
    2296 !  which can later be used to make the statistic files of the run:
    2297 !  if flag "callstats" from callphys.def is .true.)
    2298 
    2299 
    2300       call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
    2301       call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
    2302       call wstats(ngrid,"fluxsurf_lw",                               &
    2303                   "Thermal IR radiative flux to surface","W.m-2",2,  &
    2304                   fluxsurf_lw)
    2305       call wstats(ngrid,"fluxtop_lw",                                &
    2306                   "Thermal IR radiative flux to space","W.m-2",2,    &
    2307                   fluxtop_lw)
    2308 
    2309       ! call wstats(ngrid,"fluxsurf_sw",                               &
    2310       !             "Solar radiative flux to surface","W.m-2",2,       &
    2311       !             fluxsurf_sw_tot)
    2312       ! call wstats(ngrid,"fluxtop_sw",                                &
    2313       !             "Solar radiative flux to space","W.m-2",2,         &
    2314       !             fluxtop_sw_tot)
    2315 
    2316       call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
    2317       call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
    2318       call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
    2319       ! call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
    2320       ! call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
    2321       call wstats(ngrid,"p","Pressure","Pa",3,pplay)
    2322       call wstats(ngrid,"emis","Emissivity","",2,emis)
    2323       call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
    2324       call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
    2325       call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
    2326       call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
    2327       call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
    2328 
    2329       if (tracer) then
    2330          do iq=1,nq
    2331             call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
    2332             call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',&
    2333                         'kg m^-2',2,qsurf(1,iq) )
    2334             call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',  &
    2335                         'kg m^-2',2,qcol(1,iq) )
    2336             ! call wstats(ngrid,trim(noms(iq))//'_reff',trim(noms(iq))//'_reff',&
    2337             !             'm',3,reffrad(1,1,iq)) ! bug (2): Subscript #3 of the array REFFRAD has value 2 which is greater than the upper bound of 1
    2338          end do
    2339 
    2340       endif ! end of 'tracer'
    2341 
    2342       if(lastcall.and.callstats) then
    2343          write (*,*) "Writing stats..."
    2344          call mkstats(ierr)
    2345       endif
    2346 
    2347 #ifndef MESOSCALE
    23482187
    23492188!------------------------------------------------------------------------------
    … …  
    23572196!------------------------------------------------------------------------------
    23582197
     2198      !-------- General 1D variables
     2199
    23592200      call write_output("Ls","solar longitude","deg",zls*180./pi)
    23602201      ! call write_output("Lss","sub solar longitude","deg",zlss*180./pi)
    … …  
    23622203      call write_output("Declin","solar declination","deg",declin*180./pi)
    23632204      call write_output("dist_star","dist_star","AU",dist_star)
     2205      call write_output("globave","surf press","Pa",globave)
     2206
     2207      !-------- General 2D variables
    23642208
    23652209      call write_output("tsurf","Surface temperature","K",tsurf)
    23662210      call write_output("ps","Surface pressure","Pa",ps)
    23672211      call write_output("emis","Emissivity","",emis)
     2212      !if (grid_type == regular_lonlat) then
     2213      !     call write_output("area","Mesh area","m2", &
     2214      !                      cell_area_for_lonlat_outputs)
     2215      !   else ! unstructured grid (e.g. dynamico)
     2216      !     call write_output("area","Mesh area","m2",cell_area)
     2217      !endif
    23682218
    23692219      if (fast) then
    2370          call write_output("globave","surf press","Pa",globave)
    23712220         call write_output("fluxrad","fluxrad","W m-2",fluxrad)
    23722221         call write_output("fluxgrd","fluxgrd","W m-2",fluxgrd)
     2222         ! call write_output("dplanck","dplanck","W.s m-2 K-1",dplanck)
     2223         ! "soil" variables
    23732224         call write_output("capcal","capcal","W.s m-2 K-1",capcal)
    2374          ! call write_output("dplanck","dplanck","W.s m-2 K-1",dplanck)
    23752225         call write_output("tsoil","tsoil","K",tsoil)
    2376       else
     2226      endif
     2227
     2228      ! Total energy balance diagnostics
     2229      if(callrad)then
     2230         call write_output("ALB","Surface albedo"," ",albedo_equivalent)
     2231         call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw)
     2232         call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn)
     2233         call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw)
     2234         call write_output("GND","heat flux from ground","W m-2",fluxgrd)
     2235         if (.not.fast) then
     2236           call write_output("DYN","dynamical heat input","W m-2",fluxdyn)
     2237         endif
     2238      endif ! end of 'callrad'
     2239
     2240      !-------- General 3D variables
     2241
     2242      if (.not.fast) then
    23772243         if (check_physics_outputs) then
    23782244            ! Check the validity of updated fields at the end of the physics step
    … …  
    23932259      endif
    23942260
    2395       ! Total energy balance diagnostics
    2396       if(callrad)then
    2397          call write_output("ALB","Surface albedo"," ",albedo_equivalent)
    2398          call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw)
    2399          call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn)
    2400          call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw)
    2401          call write_output("GND","heat flux from ground","W m-2",fluxgrd)
    2402          call write_output("DYN","dynamical heat input","W m-2",fluxdyn)
    2403       endif ! end of 'callrad'
    2404 
    24052261      if(enertest) then
    24062262         if (calldifv) then
    … …  
    24212277      endif ! end of 'enertest'
    24222278
    2423          ! Diagnostics of optical thickness
    2424          ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
    2425          if (diagdtau) then
     2279      ! Diagnostics of optical thickness
     2280      ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
     2281      if (diagdtau) then
    24262282            do nw=1,L_NSPECTV
    24272283               write(str2,'(i2.2)') nw
    … …  
    24322288               call write_output('dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',int_dtaui(:,nlayer:1:-1,nw))
    24332289            enddo
    2434          endif
    2435 
    2436          ! Temporary inclusions for heating diagnostics.
    2437          call write_output("zdtsw","SW heating","T s-1",zdtsw)
    2438          call write_output("zdtlw","LW heating","T s-1",zdtlw)
    2439          call write_output("dtrad","radiative heating","K s-1",dtrad)
    2440          call write_output("zdtdyn","Dyn. heating","T s-1",zdtdyn)
    2441 
    2442          ! For Debugging.
    2443          !call write_output('rnat','Terrain type',' ',real(rnat))
     2290      endif
     2291
     2292      ! Temporary inclusions for heating diagnostics.
     2293      if (.not.fast) then
     2294        call write_output("zdtsw","SW heating","T s-1",zdtsw)
     2295        call write_output("zdtlw","LW heating","T s-1",zdtlw)
     2296        call write_output("dtrad","radiative heating","K s-1",dtrad)
     2297        call write_output("zdtdyn","Dyn. heating","T s-1",zdtdyn)
     2298      endif
     2299
     2300      ! For Debugging.
     2301      !call write_output('rnat','Terrain type',' ',real(rnat))
    24442302
    24452303      ! Output tracers.
    24462304      if (tracer) then
    2447          ! call write_output("zdtc","tendancy T cond N2","K",zdtc)
    24482305
    24492306         do iq=1,nq
    2450             call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq))
    2451             ! call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    2452          !                    'kg m^-2',qsurf_hist(1,iq) )
     2307            if (.not.fast) then
     2308              call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq))
     2309            endif
    24532310            call write_output(trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
    24542311                           'kg m^-2',qcol(:,iq) )
    … …  
    24572314         enddo ! end of 'nq' loop
    24582315
    2459          !Pluto specific
    2460          call write_output('n2_iceflux','n2_iceflux',"kg m^-2 s^-1",flusurf(1,igcm_n2) )
    2461          if (haze_radproffix)then
    2462             call write_output('haze_reff','haze_reff','m',reffrad(1,1,1))
    2463          end if
     2316         ! N2 cycle
     2317         call write_output('n2_iceflux','n2_iceflux',"kg m^-2 s^-1",flusurf(:,igcm_n2) )
     2318         if (.not.fast) then
     2319            call write_output("zdtc","tendancy T cond N2","K",zdtc)
     2320         endif
     2321
     2322         ! CH4 cycle
    24642323         if (methane) then
    2465             ! call write_output("rice_ch4","ch4 ice mass mean radius","m",rice_ch4)
    2466             ! call write_output("zq1temp_ch4"," "," ",zq1temp_ch4)
    2467             ! call write_output("qsat_ch4"," "," ",qsat_ch4)
    2468             ! call write_output("qsat_ch4_l1"," "," ",qsat_ch4_l1)
    24692324
    24702325            call write_output('ch4_iceflux','ch4_iceflux',&
    2471                               "kg m^-2 s^-1",flusurf(1,igcm_ch4_ice) )
     2326                              "kg m^-2 s^-1",flusurf(:,igcm_ch4_ice) )
    24722327            call write_output("vmr_ch4","vmr_ch4","%",vmr_ch4)
     2328
    24732329            if (.not.fast) then
    24742330               call write_output("zrho_ch4","zrho_ch4","kg.m-3",zrho_ch4(:,:))
     2331               !call write_output("rice_ch4","ch4 ice mass mean radius","m",rice_ch4)
     2332               !call write_output("zq1temp_ch4"," "," ",zq1temp_ch4)
     2333               !call write_output("qsat_ch4"," "," ",qsat_ch4)
     2334               !call write_output("qsat_ch4_l1"," "," ",qsat_ch4_l1)
     2335
     2336               ! 3D Tendancies
     2337               call write_output("zdqcn2_ch4","zdq condn2 ch4","",&
     2338                           zdqc(:,:,igcm_ch4_gas))
     2339               call write_output("zdqdif_ch4","zdqdif ch4","",&
     2340                           zdqdif(:,:,igcm_ch4_gas))
     2341               call write_output("zdqsdif_ch4_ice","zdqsdif ch4","",&
     2342                           zdqsdif(:,igcm_ch4_ice))
     2343               call write_output("zdqadj_ch4","zdqadj ch4","",&
     2344                           zdqadj(:,:,igcm_ch4_gas))
    24752345            endif
    24762346
    2477             ! Tendancies
    2478             call write_output("zdqcn2_ch4","zdq condn2 ch4","",&
    2479                            zdqc(:,:,igcm_ch4_gas))
    2480             call write_output("zdqdif_ch4","zdqdif ch4","",&
    2481                            zdqdif(:,:,igcm_ch4_gas))
    2482             call write_output("zdqsdif_ch4_ice","zdqsdif ch4","",&
    2483                            zdqsdif(:,igcm_ch4_ice))
    2484             call write_output("zdqadj_ch4","zdqadj ch4","",&
    2485                            zdqadj(:,:,igcm_ch4_gas))
    24862347            if (sedimentation) then
    24872348               call write_output("zdqsed_ch4","zdqsed ch4","",&
    … …  
    24902351                              zdqssed(:,igcm_ch4_gas))
    24912352            endif
     2353
    24922354            if (metcloud.and.(.not.fast)) then
    24932355               call write_output("zdtch4cloud","ch4 cloud","T s-1",&
    24942356                           zdtch4cloud)
    24952357               call write_output("zdqch4cloud","ch4 cloud","T s-1",&
    2496                            zdqch4cloud(1,1,igcm_ch4_gas))
     2358                           zdqch4cloud(:,:,igcm_ch4_gas))
    24972359            endif
    24982360
    24992361         endif
    25002362
     2363         ! CO cycle
    25012364         if (carbox) then
    25022365            ! call write_output("zdtcocloud","tendancy T cocloud","K",zdtcocloud)
    25032366            call write_output('co_iceflux','co_iceflux',&
    2504                                "kg m^-2 s^-1",flusurf(1,igcm_co_ice) )
     2367                               "kg m^-2 s^-1",flusurf(:,igcm_co_ice) )
    25052368            call write_output("vmr_co","vmr_co","%",vmr_co)
    25062369            if (.not.fast) THEN
    … …  
    25092372         endif
    25102373
     2374         ! Haze
    25112375         if (haze) then
    2512    !         call write_output("zrho_haze","zrho_haze","kg.m-3",zrho_haze(:,:))
     2376
     2377            if (haze_radproffix)then
     2378               call write_output('haze_reff','haze_reff','m',reffrad(:,:,1))
     2379            end if
     2380            !call write_output("zrho_haze","zrho_haze","kg.m-3",zrho_haze(:,:))
     2381            !call write_output("zdqhaze_col","zdqhaze col","kg/m2/s",&
     2382            !                        zdqhaze_col(:))
     2383
     2384            ! 3D Tendencies
    25132385            call write_output("zdqrho_photprec","zdqrho_photprec",&
    25142386                        "kg.m-3.s-1",zdqrho_photprec(:,:))
    … …  
    25192391            call write_output("zdqhaze_prec","zdqhaze_prec","",&
    25202392                     zdqhaze(:,:,igcm_prec_haze))
     2393            call write_output("zdqphot_ch4","zdqphot_ch4","",&
     2394                                                zdqphot_ch4(:,:))
     2395            call write_output("zdqconv_prec","zdqconv_prec","",&
     2396                                                zdqconv_prec(:,:))
     2397
    25212398            if (igcm_haze.ne.0) then
    25222399               call write_output("zdqhaze_haze","zdqhaze_haze","",&
    … …  
    25272404               endif
    25282405            endif
    2529             call write_output("zdqphot_ch4","zdqphot_ch4","",&
    2530                                                 zdqphot_ch4(:,:))
    2531             call write_output("zdqconv_prec","zdqconv_prec","",&
    2532                                                 zdqconv_prec(:,:))
    2533    !         call write_output("zdqhaze_col","zdqhaze col","kg/m2/s",
    2534    !     &                   zdqhaze_col(:))
    2535          endif
    2536 
    2537          if (aerohaze) then
    2538             call write_output("tau_col",&
    2539                "Total aerosol optical depth","opacity",tau_col)
     2406
     2407           if (aerohaze) then
     2408              call write_output("tau_col",&
     2409                 "Total aerosol optical depth","opacity",tau_col)
     2410           endif
     2411
    25402412         endif
    25412413
    25422414      endif ! end of 'tracer'
    2543 
    25442415
    25452416      ! Output spectrum.
    … …  
    25502421      endif
    25512422
    2552 #else
    2553    comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
    2554    comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
    2555    comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
    2556    if (.not.calldifv) comm_LATENT_HF(:)=0.0
    2557 
    2558    if ((tracer).and.(generic_condensation)) then
    2559 
    2560       ! If you have set generic_condensation (and not water) and you have set several GCS
    2561       ! then the outputs given to WRF will be only the ones for the last generic tracer
    2562       ! (Because it is rewritten every tracer in the loop)
    2563       ! WRF can take only one moist tracer
    2564 
    2565       do iq=1,nq
    2566          call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)
    2567 
    2568          if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
    2569 
    2570             reffrad_generic_zeros_for_wrf(:,:) = 1.
    2571 
    2572             comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer)
    2573             comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !??????
    2574             comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq)
    2575             comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas)
    2576             comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice)
    2577             ! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros !
    2578             !comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros !
    2579             comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq)
    2580             comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer)
    2581             comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer)
    2582             comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
    2583 
    2584          endif
    2585       end do ! do iq=1,nq loop on tracers
    2586 
    2587    else
    2588       comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
    2589       comm_TOTCLOUDFRAC(1:ngrid)=0.
    2590       comm_SURFRAIN(1:ngrid)=0.
    2591       comm_DQVAP(1:ngrid,1:nlayer)=0.
    2592       comm_DQICE(1:ngrid,1:nlayer)=0.
    2593       ! comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
    2594       comm_REEVAP(1:ngrid)=0.
    2595       comm_DTRAIN(1:ngrid,1:nlayer)=0.
    2596       comm_DTLSC(1:ngrid,1:nlayer)=0.
    2597       comm_RH(1:ngrid,1:nlayer)=0.
    2598 
    2599    endif ! if water, if generic_condensation, else
    2600 
    2601    comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
    2602    comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)
    2603    comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid)
    2604    comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid)
    2605    comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid)
    2606    comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid)
    2607    sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ????
    2608    comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer)
    2609 #endif
    2610 
    26112423! XIOS outputs
    26122424#ifdef CPP_XIOS

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/planete_mod.F90

    r3241 r3539  
    88  REAL,SAVE :: obliquit ! Obliquity of the planet (deg)
    99!$OMP THREADPRIVATE(apoastr,periastr,year_day,peri_day,obliquit)
    10   REAL,SAVE :: nres ! tidal resonance ratio
    1110  REAL,SAVE :: z0 ! surface roughness (m)
    1211  REAL,SAVE :: lmixmin ! mixing length
    1312  REAL,SAVE :: emin_turb ! minimal energy
    14 !$OMP THREADPRIVATE(nres,z0,lmixmin,emin_turb)
     13!$OMP THREADPRIVATE(z0,lmixmin,emin_turb)
    1514  REAL,SAVE :: coefvis
    1615  REAL,SAVE :: coefir

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/surfprop.F90

    r3483 r3539  
    7070      REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
    7171      REAL,INTENT(OUT) :: emis(ngrid)
    72       REAL,INTENT(OUT) :: therm_inertia(ngrid,nsoilmx) ! therm_inertia = inertiedat
     72      REAL,INTENT(OUT) :: therm_inertia(ngrid,nsoilmx)
    7373!-----------------------------------------------------------------------
    7474!     Local variables
    … …  
    419419           ! get depth of the different layers
    420420           ! limit diurnal / seasonal
    421            if (changetid) then
     421           if (changetid.and.methane) then
    422422              if (qsurf(ig,igcm_n2)>1.e-3) then
    423423                 emin=facls*ITN2d/volcapa*(daysec/pi)**0.5
    … …  
    426426                 emin=facls*ITCH4d/volcapa*(daysec/pi)**0.5
    427427                 tid=ITCH4d
     428              else
     429                 emin=facls*ITH2Od/volcapa*(daysec/pi)**0.5
     430                 tid=ITH2Od
     431              endif
     432           else if (changetid) then
     433              if (qsurf(ig,igcm_n2)>1.e-3) then
     434                 emin=facls*ITN2d/volcapa*(daysec/pi)**0.5
     435                 tid=ITN2d
    428436              else
    429437                 emin=facls*ITH2Od/volcapa*(daysec/pi)**0.5

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/testconserv.F90

    r3412 r3539  
    1       subroutine testconserv(ngrid,nlayer,nq,igcm1,igcm2,ptimestep, &
     1      subroutine testconserv(ngrid,nlayer,nq,pq,pdq,pqs,pdqs, &
     2          igcm1,igcm2,ptimestep, &
    23          pplev,zdq,zdqs,car1,car2)
    34
    4       use comgeomfi_h         
    55      use comcstfi_mod, only: pi, g
    66      use geometry_mod, only: latitude, longitude, cell_area
    77      USE tracer_h, only: igcm_co_ice, igcm_ch4_ice
     8      USE comgeomfi_h, only: totarea,totarea_planet
    89      implicit none
    910
    … …  
    3940      INTEGER ngrid, nlayer, nq
    4041      INTEGER igcm1, igcm2
     42      REAL pq(ngrid,nlayer,nq)
     43      REAL pdq(ngrid,nlayer,nq)
    4144      REAL zdq(ngrid,nlayer,nq)
     45      REAL pqs(ngrid,nq)
     46      REAL pdqs(ngrid,nq)
    4247      REAL zdqs(ngrid,nq)
    4348      REAL ptimestep
    … …  
    4954      INTEGER l,ig,iq
    5055      REAL masse
     56      REAL pqc(ngrid,nlayer,nq)
     57      REAL pqcs(ngrid,nq)
    5158
    5259      ! OUTPUT
    5360      REAL dWtot   
    5461      REAL dWtots
     62      REAL Wtot   
     63      REAL Wtots
    5564      REAL nconsMAX
    5665      INTEGER myig
    … …  
    5867!-----------------------------------------------------------------------
    5968
    60       write(*,*) 'TB igcm=',igcm1,igcm2
    61       write(*,*) 'TB zdq1=',zdq(1,1,igcm1)
    6269      dWtot=0.0
    6370      dWtots=0.0
     71      Wtot=0.0
     72      Wtots=0.0
    6473      nconsMAX=0.0
     74      pqc=pq+pdq*ptimestep
     75      pqcs=pqs+pdqs*ptimestep
    6576      do ig = 1, ngrid
    6677         vdifcncons(ig)=0.0
    … …  
    7081             iq    = igcm1
    7182             ! sum of atmospheric mass : kg
     83             Wtot = Wtot + masse*pqc(ig,l,iq)*cell_area(ig)
    7284             dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*cell_area(ig)
    7385             ! for each column, total mass lost per sec : kg(tracer) / m2 /s
    … …  
    7890              iq    = igcm2
    7991              dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*cell_area(ig)
     92              Wtot = Wtot + masse*pqc(ig,l,iq)*cell_area(ig)
    8093              vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq)
    8194             ENDIF
    … …  
    8396         enddo
    8497
    85              iq     = igcm1
    86              dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig)
    87              vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq)
     98         iq     = igcm1
     99         dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig)
     100         Wtots = Wtots + pqcs(ig,iq)*cell_area(ig)
     101         vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq)
    88102
    89              IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN
     103         IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN
    90104              iq     = igcm2
    91105              dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig)
     106              Wtots = Wtots + pqcs(ig,iq)*cell_area(ig)
    92107              vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq)
    93              ENDIF
     108         ENDIF
    94109             
    95              ! vdifcncons is the total amount of material that appear or
    96              ! disapear per second in the routine
    97              ! it is the non conservative factor
     110        ! vdifcncons is the total amount of material that appear or
     111        ! disapear per second in the routine
     112        ! it is the non conservative factor
    98113
    99              if(vdifcncons(ig).gt.nconsMAX)then
    100                 nconsMAX=vdifcncons(ig)
    101                 myig=ig
    102              endif
     114        if(abs(vdifcncons(ig)).gt.abs(nconsMAX))then
     115              nconsMAX=vdifcncons(ig)
     116              myig=ig
     117        endif
    103118
    104119      enddo
    105       write(*,*) 'TB Dwtot=',dWtots,dWtot
    106120
    107       dWtot  = dWtot/totarea
    108       dWtots = dWtots/totarea
    109       print*,'-------------------------------------------'
     121      dWtot  = dWtot/totarea_planet
     122      dWtots = dWtots/totarea_planet
    110123      print*,'In ',car2,' atmospheric ',car1,' change=',dWtot,' kg m-2'
    111124      print*,'In ',car2,' surface ',car1,' change  =',dWtots,' kg m-2'
    … …  
    115128        print*,'--> obtained at lat/lon=',latitude(myig)*180./pi,longitude(myig)*180./pi
    116129      ENDIF
     130      print*,'--> Total Mass ',car1,' =',Wtot+Wtots,' kg m-2'
    117131      end subroutine testconserv
    118132

  • TabularUnified trunk/LMDZ.PLUTO/libf/phypluto/testconservfast.F90

    r3412 r3539  
    1       subroutine testconservfast(ngrid,nlayer,nq,ptimestep, &
    2           pplev,zdq,zdqs,car1,car2)
     1      subroutine testconservfast(ngrid,nlayer,nq,pq,pdq,pqs,pdqs, &
     2          ptimestep,pplev,zdq,zdqs,car1,car2)
    33
    4       use comgeomfi_h         
    54      use comcstfi_mod, only: pi, g
    65      use geometry_mod, only: latitude, longitude, cell_area
     6      USE comgeomfi_h, only: totarea,totarea_planet
    77      implicit none
    88
    … …  
    3737
    3838      INTEGER ngrid, nlayer, nq
     39
     40      REAL pq(ngrid)
     41      REAL pdq(ngrid)
    3942      REAL zdq(ngrid)
     43      REAL pqs(ngrid)
     44      REAL pdqs(ngrid)
    4045      REAL zdqs(ngrid)
    4146      REAL ptimestep
    … …  
    4752      INTEGER l,ig,iq
    4853      REAL masse
     54      REAL pqc(ngrid)
     55      REAL pqcs(ngrid)
    4956
    5057      ! OUTPUT
    5158      REAL dWtot   
    5259      REAL dWtots
     60      REAL Wtot   
     61      REAL Wtots
    5362      REAL nconsMAX
    5463      INTEGER myig
    … …  
    5867      dWtot=0.0
    5968      dWtots=0.0
     69      Wtot=0.0
     70      Wtots=0.0
    6071      nconsMAX=0.0
     72      pqc=pq+pdq*ptimestep
     73      pqcs=pqs+pdqs*ptimestep
    6174      do ig = 1, ngrid
    62          vdifcncons(ig)=0.0
     75        vdifcncons(ig)=0.0
    6376
    6477        ! sum of atmospheric mass : kg
    65          dWtot = dWtot + zdq(ig)*ptimestep*cell_area(ig)*pplev(ig,1)/g
     78        Wtot = Wtot + pqc(ig)*cell_area(ig)*pplev(ig,1)/g
     79        dWtot = dWtot + zdq(ig)*ptimestep*cell_area(ig)*pplev(ig,1)/g
    6680        ! for each column, total mass lost per sec : kg(tracer) / m2 /s
    67          vdifcncons(ig)=vdifcncons(ig) + pplev(ig,1)/g*zdq(ig)
     81        vdifcncons(ig)=vdifcncons(ig) + pplev(ig,1)/g*zdq(ig)
    6882 
    69          dWtots = dWtots + zdqs(ig)*ptimestep*cell_area(ig)
    70          vdifcncons(ig)=vdifcncons(ig)+zdqs(ig)
     83        dWtots = dWtots + zdqs(ig)*ptimestep*cell_area(ig)
     84        Wtots = Wtots + pqcs(ig)*cell_area(ig)
     85        vdifcncons(ig)=vdifcncons(ig)+zdqs(ig)
    7186
    7287             
    73              ! vdifcncons is the total amount of material that appear or
    74              ! disapear per second in the routine
    75              ! it is the non conservative factor
     88        ! vdifcncons is the total amount of material that appear or
     89        ! disapear per second in the routine
     90        ! it is the non conservative factor
    7691
    77          if(vdifcncons(ig).gt.nconsMAX)then
     92        if(abs(vdifcncons(ig)).gt.abs(nconsMAX))then
    7893              nconsMAX=vdifcncons(ig)
    7994              myig=ig
    80          endif
     95        endif
    8196
    8297      enddo
    8398
    84       dWtot  = dWtot/totarea
    85       dWtots = dWtots/totarea
    86       print*,'-------------------------------------------'
     99      dWtot  = dWtot/totarea_planet
     100      dWtots = dWtots/totarea_planet
    87101      print*,'In ',car2,' atmospheric ',car1,' change=',dWtot,' kg m-2'
    88102      print*,'In ',car2,' surface ',car1,' change  =',dWtots,' kg m-2'
    … …  
    92106        print*,'--> obtained at lat/lon=',latitude(myig)*180./pi,longitude(myig)*180./pi
    93107      ENDIF
     108      print*,'--> Total Mass ',car1,' =',(Wtot+Wtots)/totarea_planet,' kg m-2'
    94109      end subroutine testconservfast
    95110
Note: See TracChangeset for help on using the changeset viewer.