Changeset 883 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Feb 15, 2013, 2:55:26 PM (12 years ago)
Author:
tnavarro
Message:

corrected problems with update in pressure levels, corrected bug without scavenging

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r835 r883  
    128128#include "tracer.h"
    129129#include "nlteparams.h"
     130#include "comvert.h"
    130131
    131132#include "chimiedata.h"
    … …  
    306307
    307308      real co2col(ngridmx)        ! CO2 column
    308       REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
     309      ! pplev and pplay are dynamical inputs and must not be modified in the physics.
     310      ! instead, use zplay and zplev :
     311      REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
    309312      REAL zstress(ngrid), cd
    310313      real hco2(nqmx),tmean, zlocal(nlayermx)
    … …  
    511514      end if
    512515
     516c     Initialize pressure levels
     517c     ~~~~~~~~~~~~~~~~~~
     518      zplev(:,:) = pplev(:,:)
     519      zplay(:,:) = pplay(:,:)
     520      ps(:) = pplev(:,1)
     521
     522
    513523c     Compute geopotential at interlayers
    514524c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    526536      DO l=2,nlayer
    527537         DO ig=1,ngrid
    528             z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
    529             z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
     538            z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
     539            z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
    530540            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
    531541         ENDDO
    … …  
    539549      DO l=1,nlayer
    540550         DO ig=1,ngrid
    541             zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
     551            zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
    542552            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    543553         ENDDO
    … …  
    550560
    551561      if(photochem.or.callthermos) then
    552          call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
     562         call concentrations(zplay,pt,pdt,pq,pdq,ptimestep)
    553563      endif
    554564#endif
    … …  
    581591           IF(callnlte) then
    582592              if(nltemodel.eq.0.or.nltemodel.eq.1) then
    583                  CALL nltecool(ngrid,nlayer,nq,pplay,pt,pq,zdtnlte)
     593                 CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
    584594              else if(nltemodel.eq.2) then
    585595                co2vmr_gcm(1:ngrid,1:nlayer)=
    … …  
    596606     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
    597607                 
    598                  CALL nlte_tcool(ngrid,nlayer,pplay*9.869e-6,
     608                 CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
    599609     $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
    600610     $                ovmr_gcm,  zdtnlte )
    … …  
    609619c          Find number of layers for LTE radiation calculations
    610620           IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
    611      &          CALL nlthermeq(ngrid,nlayer,pplev,pplay)
     621     &          CALL nlthermeq(ngrid,nlayer,zplev,zplay)
    612622
    613623c          Note: Dustopacity.F has been transferred to callradite.F
    … …  
    622632
    623633           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
    624      $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
     634     $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
    625635     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
    626636     $     tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
    … …  
    635645     &             " tau(",f4.0," Pa) =",f9.6)')
    636646     &            odpref,tauref(igout),
    637      &            odpref,tau(igout,1)*odpref/pplev(igout,1)
     647     &            odpref,tau(igout,1)*odpref/zplev(igout,1)
    638648c          ---------------------------------------------------------
    639649c          Call slope parameterization for direct and scattered flux
    … …  
    680690           zdtnirco2(:,:)=0
    681691           if (callnirco2) then
    682               call nirco2abs (ngrid,nlayer,pplay,dist_sol,nq,pq,
     692              call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
    683693     .                       mu0,fract,declin, zdtnirco2)
    684694           endif
    … …  
    696706c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    697707           IF(callnlte) THEN
    698               CALL blendrad(ngrid, nlayer, pplay,
     708              CALL blendrad(ngrid, nlayer, zplay,
    699709     &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
    700710           ELSE
    … …  
    742752
    743753        CALL calldrag_noro(ngrid,nlayer,ptimestep,
    744      &                 pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
     754     &                 zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
    745755 
    746756        DO l=1,nlayer
    … …  
    799809         IF (tke_heat_flux .ne. 0.) THEN
    800810             zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)*
    801      &                 (-alog(pplay(:,1)/pplev(:,1)))
     811     &                 (-alog(zplay(:,1)/zplev(:,1)))
    802812             pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
    803813         ENDIF
    … …  
    806816         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
    807817     $        ptimestep,capcal,lwrite,
    808      $        pplay,pplev,zzlay,zzlev,z0,
     818     $        zplay,zplev,zzlay,zzlev,z0,
    809819     $        pu,pv,zh,pq,tsurf,emis,qsurf,
    810820     $        zdum1,zdum2,zdh,pdq,zflubid,
    … …  
    879889     $ zzlev,zzlay,
    880890     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
    881      $ pplay,pplev,pphi,zpopsk,
     891     $ zplay,zplev,pphi,zpopsk,
    882892     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
    883893     $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
    … …  
    931941
    932942         CALL convadj(ngrid,nlayer,nq,ptimestep,
    933      $                pplay,pplev,zpopsk,lmax_th,
     943     $                zplay,zplev,zpopsk,lmax_th,
    934944     $                pu,pv,zh,pq,
    935945     $                pdu,pdv,zdh,pdq,
    … …  
    969979         co2ice = co2ice * 10000.
    970980      ENDIF
     981     
     982     
     983      pdpsrf(:) = 0
    971984
    972985      IF (callcond) THEN
    973986         CALL newcondens(ngrid,nlayer,nq,ptimestep,
    974      $              capcal,pplay,pplev,tsurf,pt,
     987     $              capcal,zplay,zplev,tsurf,pt,
    975988     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
    976989     $              co2ice,albedo,emis,
    … …  
    9991012         ENDIF ! of IF (tracer)
    10001013
     1014        ! update surface pressure
     1015        DO ig=1,ngrid
     1016          ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
     1017        ENDDO
     1018       
     1019        ! update pressure levels
     1020        DO l=1,nlayer
     1021         DO ig=1,ngrid
     1022          zplay(ig,l) = aps(l) + bps(l)*ps(ig)
     1023          zplev(ig,l) = ap(l)  + bp(l)*ps(ig)
     1024         ENDDO
     1025        ENDDO
     1026        zplev(:,l) = 0.
     1027       
     1028        ! update layers altitude
     1029        DO l=2,nlayer
     1030          DO ig=1,ngrid
     1031           z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
     1032           z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
     1033           zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
     1034          ENDDO
     1035        ENDDO
     1036     
    10011037      ENDIF  ! of IF (callcond)
     1038     
     1039
    10021040
    10031041c-----------------------------------------------------------------------
    … …  
    10161054
    10171055           call watercloud(ngrid,nlayer,ptimestep,
    1018      &                pplev,pplay,pdpsrf,zzlay, pt,pdt,
     1056     &                zplev,zplay,pdpsrf,zzlay, pt,pdt,
    10191057     &                pq,pdq,zdqcloud,zdtcloud,
    10201058     &                nq,tau,tauscaling,rdust,rice,nuice,
    … …  
    10371075
    10381076! increment dust and ccn masses and numbers
     1077! We need to check that we have Nccn & Ndust > 0
     1078! This is due to single precision rounding problems
    10391079           if (microphys) then
     1080             pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
     1081     &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
     1082     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
     1083             pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
     1084     &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
     1085     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
     1086             where (pq(:,:,igcm_ccn_mass) +
     1087     &       ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
     1088               pdq(:,:,igcm_ccn_mass) =
     1089     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
     1090               pdq(:,:,igcm_ccn_number) =
     1091     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
     1092             end where
     1093             where (pq(:,:,igcm_ccn_number) +
     1094     &       ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
     1095               pdq(:,:,igcm_ccn_mass) =
     1096     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
     1097               pdq(:,:,igcm_ccn_number) =
     1098     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
     1099             end where
     1100           endif
     1101
     1102           if (scavenging) then
    10401103             pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
    10411104     &         pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
    … …  
    10441107     &         pdq(1:ngrid,1:nlayer,igcm_dust_number) +
    10451108     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
    1046              pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
    1047      &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
    1048      &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
    1049              pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
    1050      &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
    1051      &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
    1052            endif
    1053 
    1054 ! We need to check that we have Nccn & Ndust > 0
    1055 ! This is due to single precision rounding problems
    1056          if (scavenging) then
    1057 
    1058            where (pq(:,:,igcm_ccn_mass) +
    1059      &      ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
    1060               pdq(:,:,igcm_ccn_mass) =
    1061      &        - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
    1062               pdq(:,:,igcm_ccn_number) =
    1063      &        - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    1064            end where
    1065            where (pq(:,:,igcm_ccn_number) +
    1066      &      ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
    1067               pdq(:,:,igcm_ccn_mass) =
    1068      &        - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
    1069               pdq(:,:,igcm_ccn_number) =
    1070      &        - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    1071            end where
    1072            where (pq(:,:,igcm_dust_mass) +
    1073      &      ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
    1074               pdq(:,:,igcm_dust_mass) =
    1075      &        - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    1076               pdq(:,:,igcm_dust_number) =
    1077      &        - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    1078            end where
    1079            where (pq(:,:,igcm_dust_number) +
    1080      &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
    1081               pdq(:,:,igcm_dust_mass) =
    1082      &        - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    1083               pdq(:,:,igcm_dust_number) =
    1084      &        - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    1085            end where
    1086            
    1087          endif ! of if scavenging
     1109             where (pq(:,:,igcm_dust_mass) +
     1110     &       ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
     1111               pdq(:,:,igcm_dust_mass) =
     1112     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
     1113               pdq(:,:,igcm_dust_number) =
     1114     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1115             end where
     1116             where (pq(:,:,igcm_dust_number) +
     1117     &       ptimestep*pdq(:,:,igcm_dust_number) < 0.)
     1118               pdq(:,:,igcm_dust_mass) =
     1119     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
     1120               pdq(:,:,igcm_dust_number) =
     1121     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1122             end where
     1123           endif ! of if scavenging
    10881124                     
    10891125
    … …  
    10971133c        ----------
    10981134         IF(callddevil) then
    1099            call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
     1135           call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
    11001136     &                zdqdev,zdqsdev)
    11011137 
    … …  
    11271163
    11281164           call callsedim(ngrid,nlayer, ptimestep,
    1129      &                pplev,zzlev, zzlay, pt, rdust, rice,
     1165     &                zplev,zzlev, zzlay, pt, rdust, rice,
    11301166     &                rsedcloud,rhocloud,
    11311167     &                pq, pdq, zdqsed, zdqssed,nq,
    … …  
    11621198
    11631199!           dust and ice surface area
    1164             call surfacearea(ngrid, nlayer, ptimestep, pplay, zzlay,
     1200            call surfacearea(ngrid, nlayer, ptimestep, zplay, zzlay,
    11651201     $                       pt, pq, pdq, nq,
    11661202     $                       rdust, rice, tau, tauscaling,
    11671203     $                       surfdust, surfice)
    11681204!           call photochemistry
    1169             call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
     1205            call calchim(ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
    11701206     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
    11711207     $                   zdqcloud,zdqscloud,tauref,co2ice,
    … …  
    12331269
    12341270      if (callthermos) then
    1235         call thermosphere(pplev,pplay,dist_sol,
     1271        call thermosphere(zplev,zplay,dist_sol,
    12361272     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
    12371273     &     pt,pq,pu,pv,pdt,pdq,
    … …  
    12751311         if (caps.and.(obliquit.lt.27.)) then
    12761312           ! NB: Updated surface pressure, at grid point 'ngrid', is
    1277            !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
     1313           !     ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
    12781314           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
    1279      &                     (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
     1315     &                     (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
    12801316         endif
    12811317#endif
    … …  
    13421378            zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
    13431379          ENDDO
    1344         ENDDO
    1345       ENDDO
    1346 
    1347       ! surface pressure
    1348       DO ig=1,ngrid
    1349           ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
    1350       ENDDO
    1351 
    1352       ! pressure
    1353       DO l=1,nlayer
    1354         DO ig=1,ngrid
    1355              zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
    1356              zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
    13571380        ENDDO
    13581381      ENDDO
    … …  
    14131436         WRITE(*,'(2f10.5,2f15.5)')
    14141437     s   tsurf(igout),zdtsurf(igout)*ptimestep,
    1415      s   pplev(igout,1),pdpsrf(igout)*ptimestep
     1438     s   zplev(igout,1),pdpsrf(igout)*ptimestep
    14161439         WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
    14171440         WRITE(*,'(i4,6f10.5)') (l,
    … …  
    15451568                 mtot(ig) = mtot(ig) +
    15461569     &                      zq(ig,l,igcm_h2o_vap) *
    1547      &                      (pplev(ig,l) - pplev(ig,l+1)) / g
     1570     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
    15481571                 icetot(ig) = icetot(ig) +
    15491572     &                        zq(ig,l,igcm_h2o_ice) *
    1550      &                        (pplev(ig,l) - pplev(ig,l+1)) / g
     1573     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
    15511574c                Computing abs optical depth at 825 cm-1 in each
    15521575c                  layer to simulate NEW TES retrieval
    … …  
    15561579                 opTES(ig,l)= 0.75 * Qabsice *
    15571580     &             zq(ig,l,igcm_h2o_ice) *
    1558      &             (pplev(ig,l) - pplev(ig,l+1)) / g
     1581     &             (zplev(ig,l) - zplev(ig,l+1)) / g
    15591582     &             / (rho_ice * rice(ig,l) * (1.+nuice_ref))
    15601583                 tauTES(ig)=tauTES(ig)+ opTES(ig,l)
    … …  
    15631586c               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
    15641587             enddo
    1565              call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
     1588             call watersat(ngridmx*nlayermx,zt,zplay,zqsat)
    15661589             satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)
    15671590
    … …  
    15741597                    Nccntot(ig) = Nccntot(ig) +
    15751598     &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
    1576      &              *(pplev(ig,l) - pplev(ig,l+1)) / g
     1599     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
    15771600                    Mccntot(ig) = Mccntot(ig) +
    15781601     &              zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
    1579      &              *(pplev(ig,l) - pplev(ig,l+1)) / g
     1602     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
    15801603cccc Column integrated effective ice radius
    15811604cccc is weighted by total ice surface area (BETTER than total ice mass)
    … …  
    15831606     &                      tauscaling(ig) *
    15841607     &                      zq(ig,l,igcm_ccn_number) *
    1585      &                      (pplev(ig,l) - pplev(ig,l+1)) / g *
     1608     &                      (zplev(ig,l) - zplev(ig,l+1)) / g *
    15861609     &                      rice(ig,l) * rice(ig,l)*  (1.+nuice_ref)
    15871610                 enddo
    … …  
    15961619                 rave(ig) = rave(ig) +
    15971620     &                      zq(ig,l,igcm_h2o_ice) *
    1598      &                      (pplev(ig,l) - pplev(ig,l+1)) / g *
     1621     &                      (zplev(ig,l) - zplev(ig,l+1)) / g *
    15991622     &                      rice(ig,l) * (1.+nuice_ref)
    16001623                 enddo
    … …  
    16411664     &                "m.s-1",3,pw)
    16421665        call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
    1643         call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
     1666        call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
    16441667          call wstats(ngrid,"q2",
    16451668     &                "Boundary layer eddy kinetic energy",
    … …  
    17951818                      do ig=1,ngrid                                 
    17961819                         colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
    1797      $                                  *(pplev(ig,l)-pplev(ig,l+1))
     1820     $                                  *(zplev(ig,l)-zplev(ig,l+1))
    17981821     $                                  *6.022e22/(mmol(iq)*g)       
    17991822                      end do                                       
    … …  
    18481871        dustot(ig) = dustot(ig) +
    18491872     &               zq(ig,l,igcm_dust_mass)
    1850      &               *  (pplev(ig,l) - pplev(ig,l+1)) / g
     1873     &               *  (zplev(ig,l) - zplev(ig,l+1)) / g
    18511874       enddo
    18521875      enddo
    … …  
    19451968           do ig=1,ngrid
    19461969             co2col(ig)=co2col(ig)+
    1947      &                  zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
     1970     &                  zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
    19481971           enddo
    19491972         enddo
    … …  
    22172240           do ig=1,ngrid
    22182241             co2col(ig)=co2col(ig)+
    2219      &                  zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
     2242     &                  zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
    22202243         enddo
    22212244         enddo
    … …  
    22852308               opTES(1,l)= 0.75 * Qabsice *
    22862309     &             zq(1,l,igcm_h2o_ice) *
    2287      &             (pplev(1,l) - pplev(1,l+1)) / g
     2310     &             (zplev(1,l) - zplev(1,l+1)) / g
    22882311     &             / (rho_ice * rice(1,l) * (1.+nuice_ref))
    22892312               tauTES=tauTES+ opTES(1,l)
    … …  
    23052328           do l=1,nlayer
    23062329             mtot = mtot +  zq(1,l,igcm_h2o_vap)
    2307      &                 * (pplev(1,l) - pplev(1,l+1)) / g
     2330     &                 * (zplev(1,l) - zplev(1,l+1)) / g
    23082331             icetot = icetot +  zq(1,l,igcm_h2o_ice)
    2309      &                 * (pplev(1,l) - pplev(1,l+1)) / g
     2332     &                 * (zplev(1,l) - zplev(1,l+1)) / g
    23102333           end do
    23112334           h2otot = h2otot+mtot+icetot
    … …  
    23292352             rave = rave + tauscaling(1) *
    23302353     &              zq(1,l,igcm_ccn_number) *
    2331      &              (pplev(1,l) - pplev(1,l+1)) / g *
     2354     &              (zplev(1,l) - zplev(1,l+1)) / g *
    23322355     &              rice(1,l) * rice(1,l)*  (1.+nuice_ref)
    23332356             enddo
    … …  
    23352358
    23362359              Nccntot= 0
    2337               call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
     2360              call watersat(ngridmx*nlayermx,zt,zplay,zqsat)
    23382361               do l=1,nlayermx
    23392362                   Nccntot = Nccntot +
    23402363     &              zq(1,l,igcm_ccn_number)*tauscaling(1)
    2341      &              *(pplev(1,l) - pplev(1,l+1)) / g
     2364     &              *(zplev(1,l) - zplev(1,l+1)) / g
    23422365                   satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
    23432366                   satu(1,l) = (max(satu(1,l),float(1))-1)
    23442367!     &                      * zq(1,l,igcm_h2o_vap) *
    2345 !     &                      (pplev(1,l) - pplev(1,l+1)) / g
     2368!     &                      (zplev(1,l) - zplev(1,l+1)) / g
    23462369               enddo
    23472370               call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
    … …  
    23642387             do l=1,nlayer
    23652388               rave = rave + zq(1,l,igcm_h2o_ice)
    2366      &              * (pplev(1,l) - pplev(1,l+1)) / g
     2389     &              * (zplev(1,l) - zplev(1,l+1)) / g
    23672390     &              *  rice(1,l) * (1.+nuice_ref)
    23682391             enddo
    … …  
    23952418
    23962419
    2397          zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
     2420         zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
    23982421
    23992422         do l=2,nlayer-1
    … …  
    24022425     &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
    24032426              zlocal(l)= zlocal(l-1)
    2404      &        -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g
     2427     &        -log(zplay(1,l)/zplay(1,l-1))*rnew(1,l)*tmean/g
    24052428         enddo
    24062429         zlocal(nlayer)= zlocal(nlayer-1)-
    2407      &                   log(pplay(1,nlayer)/pplay(1,nlayer-1))*
     2430     &                   log(zplay(1,nlayer)/zplay(1,nlayer-1))*
    24082431     &                   rnew(1,nlayer)*tmean/g
    24092432
Note: See TracChangeset for help on using the changeset viewer.