Changeset 3237 for trunk/LMDZ.PLUTO.old


Ignore:
Timestamp:
Feb 27, 2024, 12:04:22 PM (9 months ago)
Author:
afalco
Message:

Pluto PCM:
Import methane and ch4 clouds from pluto.old
AF

Location:
trunk/LMDZ.PLUTO.old/libf/phypluto
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO.old/libf/phypluto/physiq.F

    r3175 r3237  
    66     *            pw,
    77     *            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
    8      
     8
    99      use radinc_h, only : naerkind
    10       use planet_h         
     10      use planet_h
    1111      USE ioipsl_getincom
    1212      use comgeomfi_h
    … …  
    4242!           - "startfi", "histfi" if it's time
    4343!           - Saving statistics if "callstats = .true."
    44 !           - Output any needed variables in "diagfi" 
     44!           - Output any needed variables in "diagfi"
    4545!     10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
    4646!
    … …  
    6868!    pt(ngrid,nlayer)      Temperature (K)
    6969!    pq(ngrid,nlayer,nq)   Advected fields
    70 !    pudyn(ngrid,nlayer)    \ 
     70!    pudyn(ngrid,nlayer)    \
    7171!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
    7272!    ptdyn(ngrid,nlayer)     / corresponding variables
    … …  
    8888!           Frederic Hourdin    15/10/93
    8989!           Francois Forget     1994
    90 !           Christophe Hourdin  02/1997 
     90!           Christophe Hourdin  02/1997
    9191!           Subroutine completly rewritten by F.Forget (01/2000)
    9292!           Water ice clouds: Franck Montmessin (update 06/2003)
    9393!           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
    9494!           New correlated-k radiative scheme: R. Wordsworth (2009)
    95 !           Many specifically Martian subroutines removed: R. Wordsworth (2009)       
    96 !           Pluto version improved :T. Bertrand (2015,2020)       
     95!           Many specifically Martian subroutines removed: R. Wordsworth (2009)
     96!           Pluto version improved :T. Bertrand (2015,2020)
    9797!==================================================================
    9898
    … …  
    126126      LOGICAL firstcall,lastcall
    127127      REAL pday
    128       REAL ptime 
     128      REAL ptime
    129129      logical tracerdyn
    130130
    … …  
    139139c Local saved variables:
    140140c ----------------------
    141       INTEGER day_ini  ! Initial date of the run (sol since Ls=0) 
     141      INTEGER day_ini  ! Initial date of the run (sol since Ls=0)
    142142      INTEGER icount     ! counter of calls to physiq during the run.
    143143      REAL tsurf(ngridmx)            ! Surface temperature (K)
    144144      REAL tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
    145       REAL,SAVE :: tidat(ngridmx,nsoilmx)    ! thermal inertia soil 
     145      REAL,SAVE :: tidat(ngridmx,nsoilmx)    ! thermal inertia soil
    146146      REAL tidat_out(ngridmx,nlayermx)  ! thermal inertia soil but output on vertical grid
    147147      REAL tsub(ngridmx)             ! temperature of the deepest layer
    … …  
    156156      REAL,SAVE :: fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
    157157      REAL,SAVE :: qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
    158       REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy 
    159       REAL qsurf1(ngridmx,nqmx)      ! saving qsurf to calculate flux over long timescales kg.m-2 
     158      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
     159      REAL qsurf1(ngridmx,nqmx)      ! saving qsurf to calculate flux over long timescales kg.m-2
    160160      REAL flusurf(ngridmx,nqmx)     ! flux cond/sub kg.m-2.s-1
    161161      REAL,SAVE :: flusurfold(ngridmx,nqmx)  ! old flux cond/sub kg.m-2.s-1
    … …  
    165165      REAL globavenewice(nqmx)          ! globalaverage 2D ice
    166166
    167       REAL,SAVE :: ptime0    ! store the first time 
     167      REAL,SAVE :: ptime0    ! store the first time
    168168
    169169      REAL dstep
    … …  
    175175
    176176
    177       REAL stephan   
     177      REAL stephan
    178178      DATA stephan/5.67e-08/  ! Stephan Boltzman constant
    179179      SAVE stephan
    … …  
    182182      SAVE acond,bcond
    183183      REAL tcond1p4Pa
    184       DATA tcond1p4Pa/38/ 
     184      DATA tcond1p4Pa/38/
    185185
    186186c Local variables :
    … …  
    191191      REAL qsurfpal(ngridmx,nqmx)           ! qsurf after a paleoclimate step : for physdem1 and restartfi
    192192      REAL phisfipal(ngridmx)               ! geopotential after a paleoclimate step : for physdem1 and restartfi
    193       REAL oblipal                   ! change of obliquity 
    194       REAL peri_daypal               ! new periday 
    195       REAL eccpal                    ! change of eccentricity 
    196       REAL tpalnew                   ! change of time 
     193      REAL oblipal                   ! change of obliquity
     194      REAL peri_daypal               ! new periday
     195      REAL eccpal                    ! change of eccentricity
     196      REAL tpalnew                   ! change of time
    197197      REAL adjustnew                 ! change in N2 ice albedo
    198198      REAL pdaypal                   ! new pday = day_ini + step
    199199      REAL zdt_tot                   ! time range corresponding to the flux of qsurfyear
    200       REAL massacc(nqmx)             ! accumulated mass 
    201       REAL masslost(nqmx)            ! accumulated mass 
    202 
    203 !     aerosol (ice or haze) extinction optical depth  at reference wavelength 
     200      REAL massacc(nqmx)             ! accumulated mass
     201      REAL masslost(nqmx)            ! accumulated mass
     202
     203!     aerosol (ice or haze) extinction optical depth  at reference wavelength
    204204!     for the "naerkind" optically active aerosols:
    205       REAL aerosol(ngridmx,nlayermx,naerkind) 
    206 
    207       CHARACTER*80 fichier 
     205      REAL aerosol(ngridmx,nlayermx,naerkind)
     206
     207      CHARACTER*80 fichier
    208208      INTEGER l,ig,ierr,igout,iq,i, tapphys
    209209      INTEGER lecttsoil     ! lecture of tsoil from proftsoil
    … …  
    218218      real,save :: fluxtop_dn(ngridmx)
    219219      real fluxabs_sw(ngridmx) ! absorbed shortwave radiation
    220       real fluxdyn(ngridmx)    ! horizontal heat transport by dynamics     
     220      real fluxdyn(ngridmx)    ! horizontal heat transport by dynamics
    221221
    222222      REAL zls                       ! solar longitude (rad)
    223       REAL zfluxuv                     ! Lyman alpha flux at 1AU 
     223      REAL zfluxuv                     ! Lyman alpha flux at 1AU
    224224                                     !  ph/cm2/s
    225225      REAL zday                      ! date (time since Ls=0, in martian days)
    226       REAL saveday                      ! saved date 
    227       SAVE saveday   
     226      REAL saveday                      ! saved date
     227      SAVE saveday
    228228      REAL savedeclin                   ! saved declin
    229229      SAVE savedeclin
    230       !REAL adjust                      ! adjustment for surface albedo 
     230      !REAL adjust                      ! adjustment for surface albedo
    231231      REAL zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
    232232      REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
    … …  
    275275      REAL zdqscofast(ngridmx)    ! used only for fast model nogcm
    276276      REAL zdqcofast(ngridmx)    ! used only for fast model nogcm
    277       REAL zdqflow(ngridmx,nqmx) 
     277      REAL zdqflow(ngridmx,nqmx)
    278278
    279279      REAL zdteuv(ngridmx,nlayermx)    ! (K/s)
    … …  
    283283      real zdqmoldiff(ngridmx,nlayermx,nqmx)
    284284
    285       ! Haze relatated tendancies 
    286       REAL zdqhaze(ngridmx,nlayermx,nqmx) 
    287       REAL zdqprodhaze(ngridmx,nqmx) 
    288       REAL zdqsprodhaze(ngridmx) 
    289       REAL zdqhaze_col(ngridmx) 
     285      ! Haze relatated tendancies
     286      REAL zdqhaze(ngridmx,nlayermx,nqmx)
     287      REAL zdqprodhaze(ngridmx,nqmx)
     288      REAL zdqsprodhaze(ngridmx)
     289      REAL zdqhaze_col(ngridmx)
    290290      REAL zdqphot_prec(ngrid,nlayer)
    291291      REAL zdqphot_ch4(ngrid,nlayer)
    292292      REAL zdqconv_prec(ngrid,nlayer)
    293       REAL zdq_source(ngridmx,nlayermx,nqmx) 
    294       ! Fast Haze relatated tendancies 
    295       REAL fluxbot(ngridmx)   
    296       REAL gradflux(ngridmx) 
     293      REAL zdq_source(ngridmx,nlayermx,nqmx)
     294      ! Fast Haze relatated tendancies
     295      REAL fluxbot(ngridmx)
     296      REAL gradflux(ngridmx)
    297297      REAL fluxlym_sol_bot(ngridmx)      ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
    298298      REAL fluxlym_ipm_bot(ngridmx)      ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
    299299      REAL flym_sol(ngridmx)      ! Incident Solar flux Lyman alpha ph.m-2.s-1
    300       REAL flym_ipm(ngridmx)      ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 
     300      REAL flym_ipm(ngridmx)      ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
    301301
    302302      real nconsMAX
    … …  
    336336      real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
    337337      real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
    338       save ztprevious 
     338      save ztprevious
    339339      real qtot1,qtot2            ! total aerosol mass
    340340      integer igmin, lmin
    … …  
    351351
    352352      real qcol(ngridmx,nqmx)
    353 !     Pluto adding variables 
     353!     Pluto adding variables
    354354      real vmr_ch4(ngridmx)  ! vmr ch4
    355355      real vmr_co(ngridmx)  ! vmr co
    … …  
    434434         ptime0=ptime
    435435         write (*,*) 'In physiq ptime0 =', ptime
    436          
     436
    437437c        Initialize albedo and orbital calculation
    438438c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    451451             phisfinew(ig)=phisfi(ig)-qsurf(ig,1)*g/1000. ! topo of bedrock below the ice
    452452           ENDDO
    453          endif 
     453         endif
    454454
    455455         if (conservn2) then
    … …  
    457457            call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
    458458         endif
    459                            
    460 c        initialize soil 
     459
     460c        initialize soil
    461461c        ~~~~~~~~~~~~~~~
    462462         IF (callsoil) THEN
    … …  
    499499
    500500      if (conservn2) then
    501          write(*,*) 'conservn2 iniloop' 
     501         write(*,*) 'conservn2 iniloop'
    502502         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
    503       endif 
    504      
    505       igout=ngrid/2+1 
     503      endif
     504
     505      igout=ngrid/2+1
    506506
    507507      zday=pday+ptime ! compute time, in sols (and fraction thereof)
    … …  
    543543c     Compute potential temperature zh
    544544      DO l=1,nlayer
    545          DO ig=1,ngrid 
     545         DO ig=1,ngrid
    546546            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
    547547            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    … …  
    555555c ajout de la conduction depuis la thermosphere
    556556      IF (callconduct) THEN
    557          
     557
    558558          call conduction (ngrid,nlayer,ptimestep,
    559559     $                  pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
    … …  
    582582             ENDDO
    583583          ENDDO
    584       ENDIF 
     584      ENDIF
    585585
    586586      IF (callmoldiff) THEN
    … …  
    601601         write(*,*) 'conservn2 thermo'
    602602         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
    603       endif 
     603      endif
    604604
    605605c-----------------------------------------------------------------------
    606606c    2. Compute radiative tendencies :
    607607c---------  --------------------------------------------------------------
    608 c     Saving qsurf to compute paleo flux condensation/sublimation         
     608c     Saving qsurf to compute paleo flux condensation/sublimation
    609609      DO iq=1, nq
    610610         DO ig=1,ngrid
    … …  
    614614             qsurf1(ig,iq)=qsurf(ig,iq)
    615615         ENDDO
    616       ENDDO 
     616      ENDDO
    617617
    618618      IF (callrad) THEN
    … …  
    634634               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
    635635     s         ztim1,ztim2,ztim3,mu0,fract)
    636          
     636
    637637           ELSE
    638638               CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
    … …  
    643643c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    644644           if (ngrid.ne.1) then
    645            
     645
    646646            !! Specific to change albedo of N2 so that Psurf
    647647            !! converges toward 1.4 Pa in "1989" seasons for Triton
    648648            !! converges toward 1.1 Pa in "2015" seasons for Pluto
    649             if (convergeps) then     
    650               if (triton) then       
     649            if (convergeps) then
     650              if (triton) then
    651651                ! 1989 declination
    652652                if (declin*180./pi.gt.-46..and.declin*180./pi.lt.-45.
    … …  
    656656                  if (globave.gt.1.5) then
    657657                        adjust=adjust+0.005
    658                   else if (globave.lt.1.3) then 
     658                  else if (globave.lt.1.3) then
    659659                        adjust=adjust-0.005
    660660                  endif
    … …  
    665665                savedeclin=declin
    666666              else
    667                 ! Pluto : 2015 declination current epoch 
     667                ! Pluto : 2015 declination current epoch
    668668                if (declin*180./pi.gt.50.and.declin*180./pi.lt.51.
    669669     &          .and.zday.gt.saveday+10000.
    … …  
    672672                  if (globave.gt.1.2) then
    673673                        adjust=adjust+0.005
    674                   else if (globave.lt.1.) then 
     674                  else if (globave.lt.1.) then
    675675                        adjust=adjust-0.005
    676676                  endif
    … …  
    698698c          Fixed zenith angle in 1D
    699699c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    700            if(ngrid.eq.1.and.globmean1d) then   
     700           if(ngrid.eq.1.and.globmean1d) then
    701701             !global mean 1D radiative equiilium
    702702             ig=1
    … …  
    719719     &             fluxtop_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,
    720720     &             firstcall,lastcall,zzlay)
    721                    
     721
    722722c            Radiative flux from the sky absorbed by the surface (W.m-2)
    723723c            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    734734                ENDDO
    735735             ENDDO
    736              
     736
    737737           else ! corrk = F
    738738
    … …  
    754754
    755755         ENDIF ! of if(mod(icount-1,iradia).eq.0) ie radiative timestep
    756    
     756
    757757!        Transformation of the radiative tendencies:
    758758!        -------------------------------------------
    … …  
    779779         write(*,*) 'conservn2 radiat'
    780780         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
    781       endif 
     781      endif
    782782
    783783!-----------------------------------------------------------------------
    … …  
    807807     &        zdudif,zdvdif,zdhdif,zdtsdif,q2,
    808808     &        zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
    809    
     809
    810810         DO l=1,nlayer
    811811            DO ig=1,ngrid
    … …  
    825825         acond=r/lw_n2
    826826
    827          if (tracer) then 
     827         if (tracer) then
    828828          DO iq=1, nq
    829829           DO l=1,nlayer
    … …  
    856856!------------------------------------------------------------------
    857857
    858       ELSE   ! case with calldifv=F   
     858      ELSE   ! case with calldifv=F
    859859         ztim1=4.*stephan*ptimestep
    860860         DO ig=1,ngrid
    … …  
    870870
    871871c        ------------------------------------------------------------------
    872 c        Methane surface sublimation and condensation in fast model (nogcm)     
     872c        Methane surface sublimation and condensation in fast model (nogcm)
    873873c        ------------------------------------------------------------------
    874874         if ((methane).and.(fast).and.condmetsurf) THEN
    … …  
    876876           call ch4surf(ngrid,nlayer,nq,ptimestep,
    877877     &     tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf,
    878      &     zdqch4fast,zdqsch4fast) 
     878     &     zdqch4fast,zdqsch4fast)
    879879
    880880           DO ig=1,ngrid
    881               dqsurf(ig,igcm_ch4_ice)= dqsurf(ig,igcm_ch4_ice) + 
     881              dqsurf(ig,igcm_ch4_ice)= dqsurf(ig,igcm_ch4_ice) +
    882882     &                    zdqsch4fast(ig)
    883               pdq(ig,1,igcm_ch4_gas)= pdq(ig,1,igcm_ch4_gas) + 
     883              pdq(ig,1,igcm_ch4_gas)= pdq(ig,1,igcm_ch4_gas) +
    884884     &                    zdqch4fast(ig)
    885885              zdtsurf(ig)=zdtsurf(ig)+lw_ch4*zdqsch4fast(ig)/capcal(ig)
    … …  
    887887         end if
    888888c        ------------------------------------------------------------------
    889 c        CO surface sublimation and condensation in fast model (nogcm)     
     889c        CO surface sublimation and condensation in fast model (nogcm)
    890890c        ------------------------------------------------------------------
    891891         if ((carbox).and.(fast).and.condcosurf) THEN
    … …  
    893893           call cosurf(ngrid,nlayer,nq,ptimestep,
    894894     &     tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf,
    895      &     zdqcofast,zdqscofast) 
     895     &     zdqcofast,zdqscofast)
    896896
    897897           DO ig=1,ngrid
    898               dqsurf(ig,igcm_co_ice)= dqsurf(ig,igcm_co_ice) + 
     898              dqsurf(ig,igcm_co_ice)= dqsurf(ig,igcm_co_ice) +
    899899     &                    zdqscofast(ig)
    900               pdq(ig,1,igcm_co_gas)= pdq(ig,1,igcm_co_gas) + 
     900              pdq(ig,1,igcm_co_gas)= pdq(ig,1,igcm_co_gas) +
    901901     &                    zdqcofast(ig)
    902902              zdtsurf(ig)=zdtsurf(ig)+lw_co*zdqscofast(ig)/capcal(ig)
    … …  
    905905
    906906      ENDIF ! of IF (calldifv)
    907      
     907
    908908      if (conservn2) then
    909909        write(*,*) 'conservn2 diff'
    910910        call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+
    911911     &                                  dqsurf(:,1)*ptimestep)
    912       endif 
     912      endif
    913913
    914914!------------------------------------------------------------------
    … …  
    952952         ENDDO
    953953
    954          if(tracer) then 
     954         if(tracer) then
    955955           DO iq=1, nq
    956956            DO l=1,nlayer
    957957              DO ig=1,ngrid
    958                  pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) 
     958                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
    959959              ENDDO
    960960            ENDDO
    … …  
    992992c     &        pdpsrf,ptimestep)
    993993
    994          if (tracer) then 
     994         if (tracer) then
    995995           zdqc(:,:,:)=0.
    996996           zdqsc(:,:)=0.
    997997           call condens_n2(ngrid,nlayer,nq,ptimestep,
    998998     &        capcal,pplay,pplev,tsurf,pt,
    999      &        pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, 
     999     &        pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
    10001000     &        qsurf(1,igcm_n2),albedo,emis,
    10011001     &        zdtc,zdtsurfc,pdpsrf,zduc,zdvc,
    … …  
    10361036            dqsurf(ig,igcm_n2)= dqsurf(ig,igcm_n2) + zdqsc(ig,igcm_n2)
    10371037         ENDDO
    1038          
     1038
    10391039      ENDIF  ! of IF (N2cond)
    10401040
    … …  
    10431043       call testconservmass(ngrid,nlayer,pplev(:,1)+
    10441044     &      pdpsrf(:)*ptimestep,qsurf(:,1)+dqsurf(:,1)*ptimestep)
    1045       endif 
    1046  
     1045      endif
     1046
    10471047!------------------------------------------------------------------
    10481048! test n2 conservation for one tendency / pplevis not updated here w pdpsrf
    … …  
    10671067
    10681068c-----------------------------------------------------------------------
    1069 c   7. Specific parameterizations for tracers 
     1069c   7. Specific parameterizations for tracers
    10701070c   -----------------------------------------
    10711071
    1072       if (tracer) then 
     1072      if (tracer) then
    10731073
    10741074c   7a. Methane, CO, and ice
    10751075c      ---------------------------------------
    1076 c      Methane ice condensation in the atmosphere 
     1076c      Methane ice condensation in the atmosphere
    10771077c      ----------------------------------------
    10781078       zdqch4cloud(:,:,:)=0.
    … …  
    11561156c         endif  ! carbox
    11571157!------------------------------------------------------------------
    1158  
    1159 c   7b. Haze particle production 
     1158
     1159c   7b. Haze particle production
    11601160c     -------------------
    11611161       IF (haze) THEN
    1162          
     1162
    11631163          zdqphot_prec(:,:)=0.
    11641164          zdqphot_ch4(:,:)=0.
    … …  
    11701170     &                                               /ptimestep
    11711171          else
    1172            call hazecloud(ngrid,nlayer,nq,ptimestep, 
     1172           call hazecloud(ngrid,nlayer,nq,ptimestep,
    11731173     &         pplay,pplev,pq,pdq,dist_sol,mu0,zfluxuv,zdqhaze,
    11741174     &         zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
    … …  
    12061206c     -------------------
    12071207
    1208 c      ------------- 
     1208c      -------------
    12091209c      Sedimentation
    1210 c      ------------- 
    1211        IF (sedimentation) THEN 
     1210c      -------------
     1211       IF (sedimentation) THEN
    12121212         call zerophys(ngrid*nlayer*nq, zdqsed)
    12131213         call zerophys(ngrid*nq, zdqssed)
    … …  
    12171217     &        pq, pdq, zdqsed, zdqssed,nq,pphi)
    12181218
    1219          DO iq=1, nq 
     1219         DO iq=1, nq
    12201220           DO l=1,nlayer
    12211221             DO ig=1,ngrid
    … …  
    12511251!        call massn2(ngrid,nlayer,pplev,qsurf(:,1),dqsurf(:,igcm_n2),
    12521252!     &        pdpsrf,ptimestep)
    1253            
     1253
    12541254c       ---------------------------------
    12551255c       Updating tracer budget on surface (before spread of glacier)
    … …  
    12581258          DO ig=1,ngrid
    12591259            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
    1260           ENDDO 
    1261         ENDDO   
    1262 
    1263       endif !  of if (tracer) 
     1260          ENDDO
     1261        ENDDO
     1262
     1263      endif !  of if (tracer)
    12641264
    12651265      if (conservn2) then
    … …  
    12671267         call testconservmass(ngrid,nlayer,pplev(:,1)+
    12681268     &       pdpsrf(:)*ptimestep,qsurf(:,1))
    1269       endif 
    1270 
    1271       DO ig=1,ngrid 
     1269      endif
     1270
     1271      DO ig=1,ngrid
    12721272        flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)-
    12731273     &                              qsurf1(ig,igcm_n2))/ptimestep
    … …  
    12851285      ENDDO
    12861286
    1287       !! Special source of haze particle ! 
     1287      !! Special source of haze particle !
    12881288      ! todo: should be placed in haze and use tendency of n2 instead of flusurf
    12891289      IF (source_haze) THEN
    1290              call hazesource(ngrid,nlayer,nq,ptimestep, 
     1290             call hazesource(ngrid,nlayer,nq,ptimestep,
    12911291     &                       pplev,flusurf,mu0,zdq_source)
    1292              
    1293              DO iq=1, nq 
     1292
     1293             DO iq=1, nq
    12941294               DO l=1,nlayer
    12951295                 DO ig=1,ngrid
    … …  
    13131313
    13141314      ENDDO
    1315    
     1315
    13161316!   9.1 Increment Surface temperature:
    13171317!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    13181318      DO ig=1,ngrid
    1319          tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
     1319         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
    13201320      ENDDO
    1321      
     1321
    13221322      ! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
    13231323      ! Lellouch et al., 2000,2011
    … …  
    13431343         tidat_out(:,l)=tidat(:,l)
    13441344      ENDDO
    1345      
     1345
    13461346!   9.3 multiply tendencies of cond/subli for paleo loop only in the
    13471347!       last Pluto year of the simulation
    … …  
    13681368            dstep=mod(zday-day_ini-ptime0,glastep)*daysec
    13691369           else
    1370             dstep=glastep*daysec       
     1370            dstep=glastep*daysec
    13711371           endif
    13721372           zdqflow(:,:)=qsurf(:,:)
    … …  
    13831383            call testconservmass(ngrid,nlayer,pplev(:,1)+
    13841384     &       pdpsrf(:)*ptimestep,qsurf(:,1))
    1385            endif 
     1385           endif
    13861386
    13871387         endif
    … …  
    13971397      flux3 = 0
    13981398      do ig=1,ngrid
    1399             flux1 = flux1 + 
     1399            flux1 = flux1 +
    14001400     &        area(ig)*(fluxtop_dn(ig) - fluxtop_sw(ig))
    1401             flux2 = flux2 + area(ig)*fluxtop_lw(ig) 
    1402             flux3 = flux3 + area(ig)*fluxtop_dn(ig) 
     1401            flux2 = flux2 + area(ig)*fluxtop_lw(ig)
     1402            flux3 = flux3 + area(ig)*fluxtop_dn(ig)
    14031403            fluxabs_sw(ig)=fluxtop_dn(ig) - fluxtop_sw(ig)
    14041404
    … …  
    14061406      print*,'Incident solar flux, absorbed solar flux, OLR (W/m2)'
    14071407      print*, flux3/totarea, '      ',  flux1/totarea ,
    1408      & '           = ', flux2/totarea 
     1408     & '           = ', flux2/totarea
    14091409
    14101410      if(meanOLR)then
    … …  
    14211421      do ig=1,ngrid
    14221422         ts1 = ts1 + area(ig)*tsurf(ig)
    1423          ts2 = min(ts2,tsurf(ig)) 
    1424          ts3 = max(ts3,tsurf(ig)) 
     1423         ts2 = min(ts2,tsurf(ig))
     1424         ts3 = max(ts3,tsurf(ig))
    14251425      end do
    1426 !      print*,'Mean Tsurf =',ts1/totarea , ' Min Tsurf=',ts2, 
     1426!      print*,'Mean Tsurf =',ts1/totarea , ' Min Tsurf=',ts2,
    14271427!     &          ' Max Tsurf =',ts3
    1428 !      print*,'Max Tsurf=',ts3,'Min Tsurf=',ts2 
     1428!      print*,'Max Tsurf=',ts3,'Min Tsurf=',ts2
    14291429
    14301430!    -------------------------------
    … …  
    14951495
    14961496      ! pressure density
    1497       IF (.not.fast) then ! 
     1497      IF (.not.fast) then !
    14981498       DO ig=1,ngrid
    14991499        DO l=1,nlayer
    … …  
    15101510!     ---------------------------------------------------------
    15111511      call zerophys(ngrid*nq,qcol)
    1512       if(tracer.and..not.fast)then 
     1512      if(tracer.and..not.fast)then
    15131513         do iq=1,nq
    15141514            DO ig=1,ngrid
    … …  
    15271527     &                        mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
    15281528          ENDDO
    1529         ELSE 
     1529        ELSE
    15301530          DO ig=1,ngrid
    15311531!           compute vmr methane
    1532             vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)*   
     1532            vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)*
    15331533     &             g/ps(ig)*mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
    1534 !           compute density methane 
     1534!           compute density methane
    15351535            DO l=1,nlayer
    1536                zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)   
     1536               zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
    15371537            ENDDO
    1538           ENDDO 
     1538          ENDDO
    15391539        ENDIF
    15401540      endif
    … …  
    15461546     &                        mmol(igcm_n2)/mmol(igcm_co_gas)*100.
    15471547          ENDDO
    1548         ELSE 
     1548        ELSE
    15491549         DO ig=1,ngrid
    15501550!           compute vmr CO
    1551             vmr_co(ig)=qcol(ig,igcm_co_gas)*   
     1551            vmr_co(ig)=qcol(ig,igcm_co_gas)*
    15521552     &             g/ps(ig)*mmol(igcm_n2)/mmol(igcm_co_gas)*100.
    15531553!           compute density CO
    15541554            DO l=1,nlayer
    1555                zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)   
     1555               zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
    15561556            ENDDO
    1557          ENDDO 
     1557         ENDDO
    15581558        ENDIF
    15591559      endif
    … …  
    15641564       DO ig=1,ngrid
    15651565        DO l=1,nlayer
    1566                zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)   
    1567                zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)   
     1566               zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
     1567               zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
    15681568            ENDDO
    1569         ENDDO 
     1569        ENDDO
    15701570      ENDIF
    15711571
    … …  
    15981598         CLOSE(13)
    15991599      ENDIF
    1600      
     1600
    16011601
    16021602      IF (ngrid.NE.1) THEN
    … …  
    16511651                  ENDIF
    16521652                ENDDO
    1653              ENDDO                 
     1653             ENDDO
    16541654             ! Finally ensure conservation of qsurf
    16551655             DO iq=1,nq
    … …  
    17211721!        "stats")          only possible in 3D runs !
    17221722
    1723          
     1723
    17241724         IF (callstats) THEN
    1725              write (*,*) "stats have been turned off in the code. 
     1725             write (*,*) "stats have been turned off in the code.
    17261726     &                You can manage your own output in physiq.F "
    17271727             stop
    … …  
    17551755
    17561756!        ---------------------------------------------------------------------
    1757 !        3D OUTPUTS 
     1757!        3D OUTPUTS
    17581758!        ----------------------------------------------------------------------
    17591759!        output in netcdf file "DIAGFI", containing any variable for diagnostic
    17601760!        (output with  period "ecritphy", set in "run.def")
    17611761!        ----------------------------------------------------------------------
    1762        
     1762
    17631763!          if(MOD(countG3D,saveG3D).eq.0)then
    17641764!          print*,'countG3D',countG3D
    … …  
    17861786     &                  albedo)
    17871787        call WRITEDIAGFI(ngrid,"emissivite","emissivite","sans dim",2,
    1788      &                  emis) 
     1788     &                  emis)
    17891789        call WRITEDIAGFI(ngrid,"fluxtop_dn","fluxtop_dn","sans dim",2,
    1790      &                  fluxtop_dn) 
     1790     &                  fluxtop_dn)
    17911791        call WRITEDIAGFI(ngrid,"ISR","incoming stellar rad.","W m-2",
    17921792     &       2,fluxtop_dn)
    … …  
    18761876        call WRITEDIAGFI(ngrid,"zdtdif","zdtdif","K",3,zdtdif)
    18771877        call WRITEDIAGFI(ngrid,"zdtconduc","tendancy T conduc",
    1878      &                                                "K",3,zdtconduc) 
     1878     &                                                "K",3,zdtconduc)
    18791879        call WRITEDIAGFI(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
    18801880        call WRITEDIAGFI(ngrid,"zdtrad","rad heating","T s-1",3,dtrad)
    … …  
    19251925     &                    trim(noms(iq))//'_col',
    19261926     &                    'kg m^-2',2,qcol(1,iq) )
    1927  
     1927
    19281928!         call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_aero',
    19291929!     &                    trim(noms(iq))//'_aero',
    … …  
    19321932
    19331933        enddo
    1934          
     1934
    19351935         call WRITEDIAGFI(ngridmx,'haze_reff',
    19361936     &                    'haze_reff',
    … …  
    20152015
    20162016        if (haze) then
    2017        
     2017
    20182018!         call WRITEDIAGFI(ngrid,"zrho_haze","zrho_haze","kg.m-3",
    20192019!     &                                            3,zrho_haze(:,:))
    … …  
    20552055
    20562056!      ----------------------------------------------------------------------
    2057 !      1D OUTPUTS 
     2057!      1D OUTPUTS
    20582058!      ----------------------------------------------------------------------
    20592059      ELSE         ! if(ngrid.eq.1)
    20602060
    20612061       if(countG1D.eq.saveG1D)then
    2062          
     2062
    20632063!      BASIC 1D ONLY
    20642064
    … …  
    20972097          call WRITEDIAGFI(ngrid,"tidat","tidat","SI",1,tidat_out(1,:))
    20982098
    2099 !       OUTPUT OF TENDANCIES         
     2099!       OUTPUT OF TENDANCIES
    21002100!          call WRITEDIAGFI(ngrid,"zdqcloud","ch4 cloud","T s-1",
    21012101!     &                  3,zdqcloud(1,1,igcm_ch4_gas))
    … …  
    21402140        ! 1D Haze
    21412141        if (haze) then
    2142        
     2142
    21432143!         call WRITEDIAGFI(ngrid,"zrho_haze","zrho_haze","kg.m-3",
    21442144!     &                                            1,zrho_haze(:,:))
    … …  
    21732173     &         "aerosol optical depth","[]",3,aerosol(1,1,1))
    21742174        endif
    2175        
     2175
    21762176!      TRACERS 1D ONLY
    21772177        if (tracer) then

  • trunk/LMDZ.PLUTO.old/libf/phypluto/vdifc.F

    r3175 r3237  
    6767
    6868c    Traceurs :
    69       integer nq 
     69      integer nq
    7070      REAL pqsurf(ngrid,nq)
    71       real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 
    72       real pdqdif(ngrid,nlay,nq) 
    73       real pdqdifeddy(ngrid,nlay,nq) 
     71      real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
     72      real pdqdif(ngrid,nlay,nq)
     73      real pdqdifeddy(ngrid,nlay,nq)
    7474      real pdqsdif(ngrid,nq),pdqsdif1(ngrid,nq)
    75      
     75
    7676c   local:
    7777c   ------
    … …  
    106106
    107107c     variable added for N2 condensation:
    108 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
     108c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    109109      REAL hh , zhcond(ngridmx,nlayermx)
    110110c      REAL latcond,tcond1p4Pa
    … …  
    116116
    117117c    Tracers :
    118 c    ~~~~~~~ 
     118c    ~~~~~~~
    119119      INTEGER iq
    120120      REAL zq(ngridmx,nlayermx,nqmx)
    … …  
    206206c     -----------------------------------
    207207
    208         if (condensn2) then 
     208        if (condensn2) then
    209209          DO ilev=1,nlay
    210210            DO ig=1,ngrid
    … …  
    216216          DO ilev=1,nlay
    217217            DO ig=1,ngrid
    218               zhcond(ig,ilev) = 0. 
     218              zhcond(ig,ilev) = 0.
    219219            END DO
    220220          END DO
    … …  
    259259        !TB16: GCM wind for flat hemisphere
    260260        IF (phisfi(ig).eq.0.) zu2=max(2.,zu2)
    261        
     261
    262262        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
    263263        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
    … …  
    265265
    266266c    ** schema de diffusion turbulente dans la couche limite
    267 c       ---------------------------------------------------- 
     267c       ----------------------------------------------------
    268268
    269269        CALL vdif_kc(ptimestep,g,pzlev,pzlay
    … …  
    322322c      --------------------------------
    323323
    324 c    ** l'equation est 
     324c    ** l'equation est
    325325c       u(t+1) =  u(t) + dt * {(du/dt)phys}(t)  +  dt * {(du/dt)difv}(t+1)
    326326c       avec
    … …  
    330330c       donc les entrees sont /zcdv/ pour la condition a la limite sol
    331331c       et /zkv/ = Ku
    332  
     332
    333333      CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
    334334      CALL multipl(ngrid,zcdv,zb0,zb)
    … …  
    363363c      --------------------------------
    364364
    365 c    ** l'equation est 
     365c    ** l'equation est
    366366c       v(t+1) =  v(t) + dt * {(dv/dt)phys}(t)  +  dt * {(dv/dt)difv}(t+1)
    367367c       avec
    … …  
    402402c      ------------------------
    403403
    404 c    ** l'equation est 
     404c    ** l'equation est
    405405c       h(t+1) =  h(t) + dt * {(dh/dt)phys}(t)  +  dt * {(dh/dt)difv}(t+1)
    406406c       avec
    … …  
    443443c       a t + \delta t,
    444444c       c'est a dire le flux radiatif a {t + \delta t}
    445 c       + le flux turbulent a {t + \delta t} 
     445c       + le flux turbulent a {t + \delta t}
    446446c            qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
    447 c            (notation K dt = /cpp*b/)       
     447c            (notation K dt = /cpp*b/)
    448448c       + le flux dans le sol a t
    449449c       + l'evolution du flux dans le sol lorsque la temperature d'interface
    … …  
    467467c        end if
    468468      ENDDO
    469      
    470 c    ** et a partir de la temperature au sol on remonte 
     469
     470c    ** et a partir de la temperature au sol on remonte
    471471c       -----------------------------------------------
    472472
    … …  
    478478      ENDDO
    479479
    480        
     480
    481481c-----------------------------------------------------------------------
    482482c   TRACERS
    … …  
    487487c     Using the wind modified by friction for lifting and  sublimation
    488488c     ----------------------------------------------------------------
    489 !     This is computed above 
     489!     This is computed above
    490490
    491491!        DO ig=1,ngrid
    … …  
    495495!        ENDDO
    496496
    497 c       Calcul flux vertical au bas de la premiere couche (cf dust on Mars) 
     497c       Calcul flux vertical au bas de la premiere couche (cf dust on Mars)
    498498c       -----------------------------------------------------------
    499         do ig=1,ngridmx 
     499        do ig=1,ngridmx
    500500          rho(ig) = zb0(ig,1) /ptimestep
    501501!          zb(ig,1) = 0.
    … …  
    503503
    504504        call zerophys(ngrid*nq,pdqsdif)
    505         pdqdif(:,:,:)=0. 
     505        pdqdif(:,:,:)=0.
    506506
    507507
    … …  
    519519!        ENDDO
    520520
    521 c        pdqdifeddy(:,:,:)=0. 
     521c        pdqdifeddy(:,:,:)=0.
    522522cc       injection a 50 km
    523523c        DO ig=1,ngrid
    … …  
    529529c        ENDDO
    530530
    531 c      Inversion pour l'implicite sur q 
     531c      Inversion pour l'implicite sur q
    532532c       --------------------------------
    533         do iq=1,nq 
    534           CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2)) 
     533        do iq=1,nq
     534          CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
    535535
    536536          if ((methane).and.(iq.eq.igcm_ch4_gas)) then
    537 c            This line is required to account for turbulent transport 
     537c            This line is required to account for turbulent transport
    538538c            from surface (e.g. ice) to mid-layer of atmosphere:
    539539             CALL multipl(ngrid,zcdv,zb0,zb(1,1))
    … …  
    548548               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
    549549               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
    550           ENDDO 
    551  
     550          ENDDO
     551
    552552          DO ilay=nlay-1,2,-1
    553553               DO ig=1,ngrid
    … …  
    559559               ENDDO
    560560          ENDDO
    561    
     561
    562562            ! special case for methane and CO ice tracer: do not include
    563563            ! ice tracer from surface (which is set when handling
    … …  
    587587            ENDDO
    588588          endif ! of if (methane.and.(iq.eq.igcm_ch4_ice))
    589    
     589
    590590c           Calculation for turbulent exchange with the surface (for ice)
    591           IF (methane.and.(iq.eq.igcm_ch4_gas)) then   
     591          IF (methane.and.(iq.eq.igcm_ch4_gas)) then
    592592
    593593            !! calcul de la valeur de q a la surface  :
    … …  
    601601            !! Prevent CH4 condensation at the surface
    602602            if (.not.condmetsurf) then
    603                qsat_ch4=qsat_ch4*1.e6 
     603               qsat_ch4=qsat_ch4*1.e6
    604604            endif
    605605
    … …  
    625625
    626626                  zq1temp_ch4(ig)=zc(ig,1)
    627                 endif   
     627                endif
    628628
    629629             zq(ig,1,igcm_ch4_gas)=zq1temp_ch4(ig)
    630630
    631 c            TB: MODIF speciale  pour  CH4 
     631c            TB: MODIF speciale  pour  CH4
    632632             pdtsrf(ig)=pdtsrf(ig)+
    633633     &        lw_ch4*pdqsdif(ig,igcm_ch4_ice)/pcapcal(ig)
    … …  
    644644            !! Prevent CO condensation at the surface
    645645            if (.not.condcosurf) then
    646                qsat_co=qsat_co*1.e6 
     646               qsat_co=qsat_co*1.e6
    647647            endif
    648648
    … …  
    654654            END DO
    655655
    656        
     656
    657657            DO ig=1,ngrid
    658658c   -------------------------------------------------------------
    659659c           TEMPORAIRE : pour initialiser CO si glacier N2
    660660c               meme si il n'y a pas de CO disponible
    661 c             if (pqsurf(ig,igcm_n2).le.10.) then 
     661c             if (pqsurf(ig,igcm_n2).le.10.) then
    662662c   -------------------------------------------------------------
    663 c 
     663c
    664664                if ((-pdqsdif(ig,igcm_co_ice)*ptimestep)
    665665     &             .gt.(pqsurf(ig,igcm_co_ice))) then
    … …  
    671671     $            (-pdqsdif(ig,igcm_co_ice)) *ptimestep) *z1(ig)
    672672                  zq1temp_co(ig)=zc(ig,1)
    673                 endif   
    674 c             endif   
     673                endif
     674c             endif
    675675
    676676               zq(ig,1,igcm_co_gas)=zq1temp_co(ig)
    … …  
    717717      ENDDO
    718718
    719       if (tracer) then 
     719      if (tracer) then
    720720        DO iq = 1, nq
    721721          DO ilev = 1, nlay
    … …  
    730730      end if
    731731
    732 c    ** diagnostique final 
     732c    ** diagnostique final
    733733c       ------------------
    734734
Note: See TracChangeset for help on using the changeset viewer.