Changeset 728


Ignore:
Timestamp:
Jul 18, 2012, 10:54:57 AM (12 years ago)
Author:
jleconte
Message:

18/07/2012 == JL

  • New water cycle scheme:
    • largescale now in F90. Robustness increased by i) including evap inside largescale ii) computing the

condensed water amount iteratively

  • same improvements in moistadj.
  • Water thermodynamical data and saturation curves centralized in module watercommn_h
    • The saturation curves used are now Tetens formula as they are analyticaly inversible (Ts(P)-> Ps(T)).

New saturation curve yields very good agreement with the former one.

  • Saturation curves are now generalized for arbitrary water amount (not just q<<1)
  • The old watersat should be removed soon.
  • The effect of water vapor on total (surface) pressure can be taken into account by setting

mass_redistrib=.true. in callphys.def (routine mass_redistribution inspired from co2_condense in martian
model but with a different scheme as many routines evaporate/condense water vapor).

  • New cloud and precipitation scheme (JL + BC):
    • The default recovery assumption for computing the total cloud fraction has been changed (total random gave too

large cloud fractions). See totalcloudfrac.F90 for details and to change this.

  • Totalcloudfraction now set the total cloud fraction to the fraction of the

optically thickest cloud and totalcloudfrac is thus called in aeropacity.

  • Only the total cloud fraction is used to compute optical depth in aeropacity (no more effective

optical depth with exponential formula).

  • 4 precipitation schemes are now available (see rain.F90 for details). The choice can be made using precip_scheme

in callphys.def. Usage of the more physically based model of Boucher et al 95 (precip_scheme=4) is recommended.
default behavior is set to the former "simple scheme" (precip_scheme=1).

  • See rain.f90 to determine the parameter to be defined in callphys.def as a function of the precipitation scheme used.
  • Physiq.F90 now written in a matricial (more F90) way.
  • Radii (H2O and CO2 cloud particles, aerosols, duts, ...) calculations now centralized in module radii_mod.F90

and work with the new aerosol scheme implemented by Laura K. Some inconsistency may remain in callsedim.


Implementation compiled with ifort and pgf90.
gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust.

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
3 added
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90

    r726 r728  
    8080      REAL tauh2so4tmp(ngridmx) ! Temporary h2so4 opacity used before scaling
    8181      ! BENJAMIN MODIFS
    82       real CLFtot,CLF1,CLF2
     82      real CLFtot
    8383      real totcloudfrac(ngridmx)
    8484      logical clearsky
    … …  
    145145!==================================================================
    146146
    147         if (iaero_co2.ne.0) then
     147      if (iaero_co2.ne.0) then
    148148           iaer=iaero_co2
    149149!       1. Initialization
    150             aerosol(:,:,iaer)=0.0
     150            aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
    151151!       2. Opacity calculation
    152152            if (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
    153                do ig=1, ngrid
    154                   do l=1,nlayer
    155                      aerosol(ig,l,iaer)=1.e-9
    156                   end do
    157                   !aerosol(ig,12,iaer)=4.0 ! single cloud layer option
    158                end do
     153               aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
     154               !aerosol(1:ngridmx,12,iaer)=4.0 ! single cloud layer option
    159155            else
    160156               DO ig=1, ngrid
    … …  
    179175               ENDDO
    180176            end if ! if fixed or varying
    181         end if ! if CO2 aerosols   
     177      end if ! if CO2 aerosols   
    182178!==================================================================
    183179!     Water ice / liquid
    184180!==================================================================
    185181
    186         if (iaero_h2o.ne.0) then
     182      if (iaero_h2o.ne.0) then
    187183           iaer=iaero_h2o
    188184!       1. Initialization
    189             aerosol(:,:,iaer)=0.0
     185            aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
    190186!       2. Opacity calculation
    191187            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
    192                do ig=1,ngrid ! to stop the rad tran bug
    193                   do l=1,nlayer
    194                      aerosol(ig,l,iaer) =1.e-9
    195                   end do
    196                end do
     188               aerosol(1:ngridmx,1:nlayermx,iaer) =1.e-9
    197189
    198190               ! put cloud at cloudlvl
    … …  
    228220
    229221            else
     222
    230223               do ig=1, ngrid
    231224                  do l=1,nlayer-1 ! to stop the rad tran bug
    232225
    233                      if(CLFvarying)then
    234                         CLF1    = max(cloudfrac(ig,l),0.01)
    235                         CLFtot  = max(totcloudfrac(ig),0.01)
    236                         CLF2    = CLF1/CLFtot
    237                         CLF2    = min(1.,CLF2)
    238                      else
    239                         CLF1=1.0
    240                         CLF2=CLFfixval
    241                      endif
    242                      
    243                      aerosol0 =                                   &
     226                     aerosol(ig,l,iaer) =                                    & !modification by BC
    244227                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
    245228                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
    246229                          !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
    247                           !( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the
    248                           !  CLF1                                   !   clear=false/pq=0 case
     230                          !( pplev(ig,l) - pplev(ig,l+1) ) / g       !   opacity in the clearsky=true and the
     231                                                                     !   clear=false/pq=0 case
    249232                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
    250                           ( pplev(ig,l) - pplev(ig,l+1) ) / g /   & 
    251                             CLF1
    252 
    253                      aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))
    254                      ! why no division in exponential?
    255                      ! because its already done in aerosol0
    256                      
    257                      aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
    258                      aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0)                     
    259                      
     233                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
     234
    260235                  enddo
    261236               enddo
    262             end if
    263         end if ! End if h2o aerosol
     237
     238               if(CLFvarying)then
     239                  call totalcloudfrac(cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngridmx,1:nlayermx,iaer))
     240                  do ig=1, ngrid
     241                     do l=1,nlayer-1 ! to stop the rad tran bug
     242                        CLFtot  = max(totcloudfrac(ig),0.01)
     243                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
     244                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
     245                     enddo
     246                  enddo
     247               else
     248                  do ig=1, ngrid
     249                     do l=1,nlayer-1 ! to stop the rad tran bug
     250                        CLFtot  = CLFfixval
     251                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
     252                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
     253                     enddo
     254                  enddo
     255              end if!(CLFvarying)
     256            endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
     257             
     258      end if ! End if h2o aerosol
    264259
    265260!==================================================================
    266261!             Dust
    267262!==================================================================
    268         if (iaero_dust.ne.0) then
     263      if (iaero_dust.ne.0) then
    269264          iaer=iaero_dust
    270 
    271265!         1. Initialization
    272             aerosol(:,:,iaer)=0.0
     266          aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
    273267         
    274             topdust=30.0 ! km  (used to be 10.0 km) LK
     268          topdust=30.0 ! km  (used to be 10.0 km) LK
    275269
    276270!       2. Opacity calculation
    … …  
    314308              ENDDO
    315309            ENDDO
    316         end if ! If dust aerosol   
     310      end if ! If dust aerosol   
    317311
    318312!==================================================================
    … …  
    321315! added by LK
    322316      if (iaero_h2so4.ne.0) then
    323           iaer=iaero_h2so4
     317         iaer=iaero_h2so4
    324318
    325319!       1. Initialization
    326             aerosol(:,:,iaer)=0.0
     320         aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
    327321
    328322
    … …  
    330324
    331325!           expfactor=0.
    332            DO l=1,nlayer-1
    333              DO ig=1,ngrid
    334 !             Typical mixing ratio profile
    335 
    336                  zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
    337                  expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
     326         DO l=1,nlayer-1
     327            DO ig=1,ngrid
     328!              Typical mixing ratio profile
     329
     330               zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
     331               expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
    338332
    339333!             Vertical scaling function
    340               aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
    341                *expfactor
    342 
    343 
    344              ENDDO
    345            ENDDO
    346             tauh2so4tmp(1:ngrid)=0.
    347               DO l=1,nlayer
    348                 DO ig=1,ngrid
    349                    tauh2so4tmp(ig) = tauh2so4tmp(ig) &
    350                           +  aerosol(ig,l,iaer)
    351                 ENDDO
    352               ENDDO
    353             DO l=1,nlayer-1
    354                DO ig=1,ngrid
    355                   aerosol(ig,l,iaer) = max(1E-20, &
     334               aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
     335
     336            ENDDO
     337         ENDDO
     338         tauh2so4tmp(1:ngrid)=0.
     339         DO l=1,nlayer
     340            DO ig=1,ngrid
     341               tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
     342            ENDDO
     343         ENDDO
     344         DO l=1,nlayer-1
     345            DO ig=1,ngrid
     346               aerosol(ig,l,iaer) = max(1E-20, &
    356347                          1 &
    357348                       *  pplev(ig,1) / preff &
    … …  
    359350                       /  tauh2so4tmp(ig))
    360351
    361               ENDDO
    362352            ENDDO
     353         ENDDO
    363354
    364355! 1/700. is assuming a "sulfurtau" of 1

  • trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90

    r726 r728  
    55!  aerosol indexes: these are initialized in iniaerosol.F and should be 0 if the
    66!                 corresponding aerosol was not activated in callphys.def
    7       integer :: iaero_co2
    8       integer :: iaero_h2o
    9       integer :: iaero_dust
    10       integer :: iaero_h2so4
     7      integer, save :: iaero_co2
     8      integer, save :: iaero_h2o
     9      integer, save :: iaero_dust
     10      integer, save :: iaero_h2so4
    1111
    1212!==================================================================

  • trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

    r726 r728  
    1414      use ioipsl_getincom
    1515      use gases_h
     16      use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad
    1617      use aerosol_mod
    1718
    … …  
    7475
    7576      REAL reffrad(ngrid,nlayer,naerkind)
    76       REAL nueffrad(ngrid,nlayer,naerkind)
     77      REAL, SAVE :: nueffrad(ngridmx,nlayermx,naerkind)
    7778
    7879!     OUTPUT
    … …  
    171172      real pqtest(ngridmx,nlayer,nq)
    172173
    173 !     are particle radii fixed?
    174       logical, save ::  radfixed
    175174      real maxrad, minrad
    176 
    177 !     Local water cloud optical properties
    178       real, parameter ::  rad_froid=35.0e-6
    179       real, parameter ::  rad_chaud=10.0e-6
    180       real, parameter ::  coef_chaud=0.13
    181       real, parameter ::  coef_froid=0.09
    182       real zfice
    183      
     175           
    184176      real CBRT
    185177      external CBRT
    … …  
    201193      qsiaer(:,:,:)=0.0
    202194      giaer(:,:,:) =0.0
    203       radfixed=.false.
    204195
    205196      if(firstcall) then
    206          if(kastprof)then
    207             radfixed=.true.
    208          endif 
    209197
    210198         call system('rm -f surf_vals_long.out')
    … …  
    212200!--------------------------------------------------
    213201!     Effective radius and variance of the aerosols
    214 
    215 
    216 !     these values will change once the microphysics gets to work
    217 !     UNLESS tracer=.false., in which case we should be working with
    218 !     a fixed aerosol layer, and be able to define reffrad in a
    219 !     .def file. To be improved!
    220 !     Generalization of aerosols added by LK
    221             do iaer = 1,naerkind
    222               if (iaer.eq.iaero_co2) then ! active CO2 aerosols 
    223                 do l=1,nlayer
    224                    do ig=1,ngrid
    225                      reffrad(ig,l,iaer)  = 1.e-4
    226                      nueffrad(ig,l,iaer) = 0.1
    227                    enddo
    228                 enddo
    229               endif
    230               if (iaer.eq.iaero_h2o) then ! active H2O aerosols
    231                 do l=1,nlayer
    232                    do ig=1,ngrid
    233                      reffrad(ig,l,iaer)  = 1.e-5
    234                      nueffrad(ig,l,iaer) = 0.1
    235                    enddo
    236                 enddo
    237               endif
    238               if(iaer.eq.iaero_dust)then ! active dust
    239                 do l=1,nlayer
    240                    do ig=1,ngrid
    241                      reffrad(ig,l,iaer)  = 1.e-5
    242                      nueffrad(ig,l,iaer) = 0.1
    243                    enddo
    244                 enddo
    245               endif
    246               if(iaer.eq.iaero_h2so4)then ! Active h2so4 aerosols
    247                 do l=1,nlayer
    248                    do ig=1,ngrid
    249                      reffrad(ig,l,iaer)  = 1.e-6
    250                      nueffrad(ig,l,iaer) = 0.1
    251                    enddo
    252                 enddo
    253               endif
    254             enddo
    255 
     202         if(naerkind.gt.4)then
     203            print*,'Code not general enough to deal with naerkind > 4 yet.'
     204            call abort
     205         endif
     206         call su_aer_radii(reffrad,nueffrad)
     207         
     208         
     209
     210!--------------------------------------------------
     211!     set up correlated k
    256212         print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
    257213         call getin("corrkdir",corrkdir)
    … …  
    305261!     Initialization on every call   
    306262
    307       do l=1,nlayer
    308          do ig=1,ngrid
    309             do iaer=1,naerkind
    310                nueffrad(ig,l,iaer) = 0.1 ! stays at 0.1
    311             enddo
    312          enddo
    313       enddo
    314 
    315 
     263!--------------------------------------------------
     264!     Effective radius and variance of the aerosols
    316265      do iaer=1,naerkind
    317266
    318        if(iaer.eq.iaero_co2)then
    319          if(radfixed)then
    320            do l=1,nlayer
    321             do ig=1,ngrid
    322                reffrad(ig,l,iaer)    = 5.e-5 ! CO2 ice
    323             enddo
    324            enddo
    325            print*,'CO2 ice particle size =',reffrad(1,1,iaer)/1.e-6,'um'
    326          else
    327            maxrad=0.0
    328            !averad=0.0
    329            minrad=1.0
    330            do l=1,nlayer
    331 
    332              !masse = (pplev(ig,l) - pplev(ig,l+1))/g
    333 
    334             do ig=1,ngrid
    335                if(tracer.and.igcm_co2_ice.gt.0)then
    336 
    337                   if(igcm_co2_ice.lt.1)then
    338                      print*,'Tracers but no CO2 ice still'
    339                      print*,'seems to be a problem...'
    340                      print*,'Aborting in callcorrk.'
    341                      stop
    342                   endif
    343                   reffrad(ig,l,iaer) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
    344                        (4*Nmix_co2*pi*rho_co2) )
    345                endif
    346                reffrad(ig,l,iaer) = max(reffrad(ig,l,iaer),1.e-6)
    347                reffrad(ig,l,iaer) = min(reffrad(ig,l,iaer),500.e-6)
    348                !averad = averad + reffrad(ig,l,iaer)*area(ig)
    349                maxrad = max(reffrad(ig,l,iaer),maxrad)
    350                minrad = min(reffrad(ig,l,iaer),minrad)
    351             enddo
    352            enddo
    353            if(igcm_co2_ice.gt.0)then
    354              print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
    355              print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
    356            endif
    357          endif
    358        endif
    359 
    360 
    361        if(iaer.eq.iaero_h2o)then
    362           if(radfixed) then
    363             do l=1,nlayer
    364                do ig=1,ngrid
    365                   !reffrad(ig,l,iaer) = 2.e-5 ! H2O ice
    366                   reffrad(ig,l,iaer) = rad_chaud ! H2O ice
    367 
    368                   zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
    369                   zfice = MIN(MAX(zfice,0.0),1.0)
    370                   reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice
    371                   nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice
    372                enddo
    373             enddo
    374             print*,'H2O ice particle size=',reffrad(1,1,iaer)/1.e-6,'um'
    375           else
    376              if(water)then
    377                maxrad=0.0
    378                minrad=1.0
    379                do l=1,nlayer
    380                 do ig=1,ngrid
    381                   reffrad(ig,l,iaer) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
    382                        (4*Nmix_h2o*pi*rho_ice) )
    383                   reffrad(ig,l,iaer) = max(reffrad(ig,l,iaer),1.e-6)
    384                   reffrad(ig,l,iaer) = min(reffrad(ig,l,iaer),100.e-6)
    385 
    386                   maxrad = max(reffrad(ig,l,iaer),maxrad)
    387                   minrad = min(reffrad(ig,l,iaer),minrad)
    388                 enddo
    389                enddo
    390                print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
    391                print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
    392              endif
    393          endif
    394        endif
    395 
     267         if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
     268            call co2_reffrad(pq,reffrad)
     269            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
     270            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
     271         end if
     272         if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
     273            call h2o_reffrad(pq,pt,reffrad,nueffrad)
     274            print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
     275            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
     276         endif
    396277         if(iaer.eq.iaero_dust)then
    397             do l=1,nlayer
    398                do ig=1,ngrid
    399                   reffrad(ig,l,iaer) = 2.e-6 ! dust
    400                enddo
    401             enddo
     278            call dust_reffrad(reffrad)
    402279            print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
    403280         endif
    404 
    405281         if(iaer.eq.iaero_h2so4)then
    406             do l=1,nlayer
    407                do ig=1,ngrid
    408                   reffrad(ig,l,iaer) = 1.e-6 ! h2so4
    409                enddo
    410             enddo
     282            call h2so4_reffrad(reffrad)
    411283            print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
    412284         endif
    413       enddo ! do iaer=1,naerkind
    414          ! added by LK
     285      end do !iaer=1,naerkind
    415286
    416287
    … …  
    641512      if(.not.kastprof)then
    642513      ! IMPORTANT: Now convert from kg/kg to mol/mol
    643       do k=1,L_LEVELS
    644          qvar(k) = qvar(k)/epsi
    645       end do
     514         do k=1,L_LEVELS
     515            qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
     516         end do
    646517      end if
    647518
    … …  
    764635            print*,'Maximum temperature is outside the radiative'
    765636            print*,'transfer kmatrix bounds, exiting.'
    766             !print*,'WARNING, OVERRIDING FOR TEST'
    767             call abort
     637            print*,'level,grid,T',k,ig,tlevrad(k)
     638            print*,'WARNING, OVERRIDING FOR TEST'
     639            !call abort
     640            tlevrad(k)=tgasmax
    768641         endif
    769642      enddo

  • trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

    r726 r728  
    99     &   , enertest                                                     &
    1010     &   , callgasvis,Continuum,H2Ocont_simple,graybody                 &
    11      &   , Nmix_co2, Nmix_h2o                                           &
    12      &   , dusttau                                                      &
    13      &   , meanOLR                                                      &
    14      &   , specOLR                                                      &
    15      &   , kastprof                                                     &
    16      &   , noradsurf                                                    &
    17      &   , Tstrat                                                       &
    18      &   , newtonian                                                    &
    19      &   , check_cpp_match                                              &
    20      &   , force_cpp                                                    &
    21      &   , tau_relax                                                    &
    22      &   , testradtimes                                                 &
     11     &   , Nmix_co2, radfixed, dusttau                                  &
     12     &   , meanOLR, specOLR                                             &
     13     &   , kastprof, noradsurf, Tstrat                                  &
     14     &   , newtonian, tau_relax, testradtimes                           &
     15     &   , check_cpp_match, force_cpp                                   &
    2316     &   , rayleigh                                                     &
    24      &   , stelbbody                                                    &
    25      &   , stelTbb                                                      &
     17     &   , stelbbody, stelTbb                                           &
    2618     &   , tplanet                                                      &
    27      &   , startype                                                     &
    28      &   , Fat1AU                                                       &
     19     &   , startype, Fat1AU                                             &
    2920     &   , nearco2cond                                                  &
    30      &   , tracer                                                       &
    31      &   , varactive                                                    &
    32      &   , varfixed                                                     &
    33      &   , satval                                                       &
     21     &   , tracer, mass_redistrib, varactive, varfixed, satval          &
    3422     &   , sedimentation,water,watercond,waterrain                      &
    35      &   , rainthreshold                                                &
    3623     &   , aeroco2,aeroh2o,aeroh2so4                                    &
    3724     &   , aerofixco2,aerofixh2o                                        &
    38      &   , hydrology                                                    &
    39      &   , sourceevol                                                   &
    40      &   , icetstep                                                     &
    41      &   , albedosnow                                                   &
    42      &   , maxicethick                                                  &
    43      &   , Tsaldiff                                                     &
    44      &   , CLFfixval                                                    &
    45      &   , CLFvarying                                                   &
     25     &   , hydrology, sourceevol, icetstep, albedosnow                  &
     26     &   , maxicethick, Tsaldiff                                        &
     27     &   , CLFfixval, CLFvarying                                        &
    4628     &   , n2mixratio                                                   &
    4729     &   , co2supsat                                                    &
    … …  
    6951      logical nearco2cond
    7052      logical tracer
     53      logical mass_redistrib
    7154      logical varactive
    7255      logical varfixed
     56      logical radfixed
    7357      logical sedimentation
    7458      logical water,watercond,waterrain
    … …  
    8771      real topdustref
    8872      real Nmix_co2
    89       real Nmix_h2o
    9073      real dusttau
    9174      real Fat1AU
    … …  
    9477      real tplanet
    9578      real satval
    96       real rainthreshold
    9779      real CLFfixval
    9880      real n2mixratio

  • trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90

    r586 r728  
    88      use radinc_h, only : naerkind
    99      use gases_h
     10      use radii_mod, only : co2_reffrad
     11      use aerosol_mod, only : iaero_co2
    1012
    1113      implicit none
    … …  
    146148      character(len=3) :: gasname
    147149
    148       real reffradmin, reffradmax
    149 
    150150      real ppco2
    151151
    … …  
    157157      call zerophys(ngridmx*nlayermx*nqmx,zq)
    158158      call zerophys(ngridmx*nlayermx,zt)
    159 
    160       !reffradmin=1.e-7
    161       !reffradmax=5.e-4
    162       !reffradmin=0.5e-7
    163       !reffradmax=0.1e-3 ! FF data
    164       reffradmin=0.1e-5
    165       reffradmax=0.1e-3 ! JB data
    166 !     improve this in the future...
    167159
    168160!     Initialisation
    … …  
    319311!     Gravitational sedimentation
    320312           
    321 !     sedimentation computed from radius computed from q
    322 !     assuming that the ice is splitted in Nmix particle
     313!     sedimentation computed from radius computed from q in module radii_mod
     314         call co2_reffrad(zq,reffrad)
     315         
    323316         do  ilay=1,nlayer
    324317            do ig=1, ngrid
    325318
    326                reff = CBRT( 3*zq(ig,ilay,i_co2ice)/( 4*Nmix_co2*pi*rho_co2 ))
    327 
    328                ! there should be a more elegant way of doing this...
    329                if(reff.lt.1.e-16) reff=1.e-16
    330 
    331                ! update reffrad for radiative transfer
    332                if(reff.lt.reffradmin)then
    333                   reffrad(ig,ilay,1)=reffradmin
    334                   !print*,'reff below optical limit'
    335                elseif(reff.gt.reffradmax)then
    336                   reffrad(ig,ilay,1)=reffradmax
    337                   !print*,'reff above optical limit'
    338                else
    339                   reffrad(ig,ilay,1)=reff
    340                endif
     319               reff = reffrad(ig,ilay,iaero_co2)
    341320
    342321               call stokes                      &

  • trunk/LMDZ.GENERIC/libf/phystd/inifis.F

    r726 r728  
    132132         write(*,*) " tracer = ",tracer
    133133
     134         write(*,*) "Run with or without atm mass update ",
     135     &      " due to tracer evaporation/condensation?"
     136         mass_redistrib=.false. ! default value
     137         call getin("mass_redistrib",mass_redistrib)
     138         write(*,*) " mass_redistrib = ",mass_redistrib
     139
    134140         write(*,*) "Diurnal cycle ?"
    135141         write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
    … …  
    377383         write(*,*)" CLFfixval = ",CLFfixval
    378384
    379          write(*,*)"Number mixing ratio of CO2 ice particles:"
    380          Nmix_co2=100000. ! default value
     385         write(*,*)"fixed radii for Cloud particles?"
     386         radfixed=.false. ! default value
     387         call getin("radfixed",radfixed)
     388         write(*,*)" radfixed = ",radfixed
     389
     390         if(kastprof)then
     391            radfixed=.true.
     392         endif 
     393
     394         write(*,*)"Number mixing ratio of CO2 ice particles:"
     395         Nmix_co2=1.e6 ! default value
    381396         call getin("Nmix_co2",Nmix_co2)
    382397         write(*,*)" Nmix_co2 = ",Nmix_co2
    383 
    384          write(*,*)"Number mixing ratio of H2O ice particles:"
    385          Nmix_h2o=10000000. ! default value
    386          call getin("Nmix_h2o",Nmix_h2o)
    387          write(*,*)" Nmix_h2o = ",Nmix_h2o
    388398
    389399!         write(*,*)"Number of radiatively active aerosols:"
    … …  
    499509         call getin("waterrain",waterrain)
    500510         write(*,*) " waterrain = ",waterrain
    501 
    502          write(*,*) "Precipitation threshold ?"
    503          rainthreshold=0.011 ! default value (Emmanuel 1997)
    504          call getin("rainthreshold",rainthreshold)
    505          write(*,*) " rainthreshold = ",rainthreshold
    506511
    507512         write(*,*) "Include surface hydrology ?"

  • trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90

    r650 r728  
    1 subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb)
    2 
    3   use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD
     1subroutine moistadj(pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
     2
     3  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water
    44
    55  implicit none
    … …  
    2424#include "comcstfi.h"
    2525
    26 !     Pre-arguments (for universal model)
    27       real pq(ngridmx,nlayermx,nqmx)       ! tracer (kg/kg)
     26      REAL pt(ngridmx,nlayermx)            ! temperature (K)
     27      REAL pq(ngridmx,nlayermx,nqmx)       ! tracer (kg/kg)
    2828      REAL pdq(ngridmx,nlayermx,nqmx)
    2929
    30       real dqmana(ngridmx,nlayermx,nqmx)   ! tendency of tracers (kg/kg.s-1)
    31       REAL dtmana(ngridmx,nlayermx)        ! temperature increment
    32 
    33 !     Arguments
    34       REAL t(ngridmx,nlayermx)       ! temperature (K)
    35       REAL q(ngridmx,nlayermx)       ! humidite specifique (kg/kg)
     30      REAL pdqmana(ngridmx,nlayermx,nqmx)  ! tendency of tracers (kg/kg.s-1)
     31      REAL pdtmana(ngridmx,nlayermx)       ! temperature increment
     32
     33!     local variables
     34      REAL zt(ngridmx,nlayermx)      ! temperature (K)
     35      REAL zq(ngridmx,nlayermx)      ! humidite specifique (kg/kg)
    3636      REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche (Pa)
    3737      REAL pplay(ngridmx,nlayermx)   ! pression au milieu de couche (Pa)
    … …  
    4949
    5050!     Local variables
    51       INTEGER i, k, iq
     51      INTEGER i, k, iq, nn
     52      INTEGER, PARAMETER :: niter = 6
    5253      INTEGER k1, k1p, k2, k2p
    5354      LOGICAL itest(ngridmx)
    … …  
    5556      REAL cp_new_t(nlayermx)
    5657      REAL cp_delta_t(nlayermx)
    57       REAL new_qb(nlayermx)
    5858      REAL v_cptj(nlayermx), v_cptjk1, v_ssig
    59       REAL v_cptt(ngridmx,nlayermx), v_p, v_t
    60       REAL v_qs(ngridmx,nlayermx), v_qsd(ngridmx,nlayermx)
     59      REAL v_cptt(ngridmx,nlayermx), v_p, v_t, zpsat
     60      REAL zqs(ngridmx,nlayermx), zdqs(ngridmx,nlayermx)
    6161      REAL zq1(ngridmx), zq2(ngridmx)
    6262      REAL gamcpdz(ngridmx,2:nlayermx)
    … …  
    9494
    9595!     GCM -----> subroutine variables
    96       DO k = 1, nlayermx
    97       DO i = 1, ngridmx
    98 
    99          q(i,k)    = pq(i,k,i_h2o)
    100 
    101          if(q(i,k).lt.0.)then
    102             q(i,k)=0.0
     96      zq(1:ngridmx,1:nlayermx)    = pq(1:ngridmx,1:nlayermx,i_h2o)+ pdq(1:ngridmx,1:nlayermx,i_h2o)*ptimestep
     97      zt(1:ngridmx,1:nlayermx)    = pt(1:ngridmx,1:nlayermx)
     98      pdqmana(1:ngridmx,1:nlayermx,1:nqmx)=0.0
     99
     100      DO k = 1, nlayermx
     101       DO i = 1, ngridmx
     102         if(zq(i,k).lt.0.)then
     103            zq(i,k)=0.0
    103104         endif
    104          DO iq = 1, nqmx
    105             dqmana(i,k,iq)=0.0
    106          ENDDO
    107       ENDDO
    108       ENDDO
    109 
    110       DO k = 1, nlayermx
    111          DO i = 1, ngridmx
    112             local_q(i,k) = q(i,k)
    113             local_t(i,k) = t(i,k)
    114             rneb(i,k) = 0.0
    115             d_ql(i,k) = 0.0
    116             d_t(i,k)  = 0.0
    117             d_q(i,k)  = 0.0
    118          ENDDO
    119          new_qb(k)=0.0
    120       ENDDO
     105       ENDDO
     106      ENDDO
     107     
     108      local_q(1:ngridmx,1:nlayermx) = zq(1:ngridmx,1:nlayermx)
     109      local_t(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
     110      rneb(1:ngridmx,1:nlayermx) = 0.0
     111      d_ql(1:ngridmx,1:nlayermx) = 0.0
     112      d_t(1:ngridmx,1:nlayermx)  = 0.0
     113      d_q(1:ngridmx,1:nlayermx)  = 0.0
    121114
    122115!     Calculate v_cptt
    … …  
    124117         DO i = 1, ngridmx
    125118            v_cptt(i,k) = RCPD * local_t(i,k)
    126             v_t = local_t(i,k)
     119            v_t = MAX(local_t(i,k),15.)
    127120            v_p = pplay(i,k)
    128121
    129             call watersat(v_t,v_p,v_qs(i,k))
    130             call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
     122            call Psat_water(v_t,v_p,zpsat,zqs(i,k))
     123            call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
    131124         ENDDO
    132125      ENDDO
    133 
    134 !     TEST: RH DIAGNOSTIC
    135 !      DO k = 1, nlayermx
    136 !         DO i = 1, ngridmx
    137 !            v_t = local_t(i,k)
    138 !            IF (v_t.LT.T_h2O_ice_liq) THEN
    139 !               print*,'RHs=',q(i,k) / v_qs(i,k)
    140 !            ELSE
    141 !               print*,'RHl=',q(i,k) / v_qs(i,k)
    142 !            ENDIF
    143 !         ENDDO
    144 !      ENDDO
    145126
    146127!     Calculate Gamma * Cp * dz: (gamma is the critical gradient)
    … …  
    149130            zdp = pplev(i,k)-pplev(i,k+1)
    150131            zdpm = pplev(i,k-1)-pplev(i,k)
    151 !         gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
    152             gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) *             &
    153                 (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
    154                 +RLVTT /(zdpm+zdp) *                            &
    155                 (v_qs(i,k-1)*zdpm + v_qs(i,k)*zdp)              &
    156                 )* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )     &
    157                 / (1.0+(v_qsd(i,k-1)*zdpm+                      &
    158                 v_qsd(i,k)*zdp)/(zdpm+zdp) )                   
     132            gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
     133                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
     134                * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
     135                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                   
    159136         ENDDO
    160137      ENDDO
    … …  
    174151      IF (k2 .GT. nlayermx) GOTO 9999
    175152      zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
    176       zsat=(local_q(i,k2-1)-v_qs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2))   &
    177          +(local_q(i,k2)-v_qs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
     153      zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2))   &
     154         +(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
    178155
    179156      IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
    … …  
    184161      IF (k2 .EQ. nlayermx) GOTO 821
    185162      k2p = k2 + 1
    186       zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
     163      zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
    187164      zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
    188165
    … …  
    194171!------------------------------------------------------ local adjustment
    195172  830 CONTINUE ! actual adjustment
     173    Do nn=1,niter
    196174      v_cptj(k1) = 0.0
    197175      zdp = pplev(i,k1)-pplev(i,k1+1)
    198       v_cptjk1 = ( (1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1))        &
    199                     + RLVTT*(local_q(i,k1)-v_qs(i,k1)) ) * zdp
    200       v_ssig = zdp * (1.0+v_qsd(i,k1))
     176      v_cptjk1 = ( (1.0+zdqs(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) + RLVTT*(local_q(i,k1)-zqs(i,k1)) ) * zdp
     177      v_ssig = zdp * (1.0+zdqs(i,k1))
    201178
    202179      k1p = k1 + 1
    … …  
    204181         zdp = pplev(i,k)-pplev(i,k+1)
    205182         v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
    206          v_cptjk1 = v_cptjk1 + zdp                                    &
    207                   * ( (1.0+v_qsd(i, k))*(v_cptt(i,k)+v_cptj(k))       &
    208                     + RLVTT*(local_q(i,k)-v_qs(i,k)) )       
    209          v_ssig = v_ssig + zdp *(1.0+v_qsd(i,k))
     183         v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))*(v_cptt(i,k)+v_cptj(k)) + RLVTT*(local_q(i,k)-zqs(i,k)) )       
     184         v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
    210185      ENDDO
    211186
    … …  
    215190         cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
    216191         cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
    217          new_qb(k) = v_qs(i,k) + v_qsd(i,k)*cp_delta_t(k)/RLVTT
    218          local_q(i,k) = new_qb(k)
     192         v_cptt(i,k)=cp_new_t(k)
     193         local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT
    219194         local_t(i,k) = cp_new_t(k) / RCPD
    220195
    221 !          print*,'v_qs in loop=',v_qs
    222 !          print*,'v_qsd in loop=',v_qsd
    223 !          print*,'new_qb in loop=',new_qb
    224 !          print*,'cp_delta_t in loop=',cp_delta_t
    225       ENDDO
     196         v_t = MAX(local_t(i,k),15.)
     197         v_p = pplay(i,k)
     198         
     199         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
     200         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
     201
     202
     203
     204!          print*,'i,k,zqs,cpdT=',i,k,zqs(i,k),cp_delta_t(k)
     205      ENDDO
     206    Enddo ! nn=1,niter
    226207
    227208
    … …  
    230211!              -- the layer about to be adjusted
    231212
    232       DO k = k1, k2
    233          v_cptt(i,k) = RCPD * local_t(i,k)
    234          v_t = local_t(i,k)
    235          v_p = pplay(i,k)
    236 
    237 !           IF (v_t.LT.t_coup) THEN
    238 !              v_qs(i,k) = qsats(v_t) / v_p
    239 !              v_qsd(i,k) = dqsats(v_t,v_qs(i,k))
    240 !           ELSE
    241 !              v_qs(i,k) = qsatl(v_t) / v_p
    242 !              v_qsd(i,k) = dqsatl(v_t,v_qs(i,k))
    243 !           ENDIF
    244 
    245          call watersat(v_t,v_p,v_qs(i,k))
    246          call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k))
    247 
    248       ENDDO
     213!      DO k = k1, k2
     214!         v_cptt(i,k) = RCPD * local_t(i,k)
     215!         v_t = local_t(i,k)
     216!         v_p = pplay(i,k)
     217!       
     218!         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
     219!         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
     220!      ENDDO
     221
    249222      DO k = 2, nlayermx
    250223         zdpm = pplev(i,k-1) - pplev(i,k)
    251224         zdp = pplev(i,k) - pplev(i,k+1)
    252 !         gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) *
    253          gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) *                       &
    254                            (v_cptt(i,k-1)*zdpm+v_cptt(i,k)*zdp)        &
    255                           +RLVTT /(zdpm+zdp) *                         &
    256                            (v_qs(i,k-1)*zdpm+v_qs(i,k)*zdp)             &
    257                          )* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )    &
    258                      / (1.0+(v_qsd(i,k-1)*zdpm+v_qsd(i,k)*zdp)         &
    259                            /(zdpm+zdp) )
     225         gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)             &
     226                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
     227                * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
     228                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                   
    260229      ENDDO
    261230
    … …  
    264233      IF (k1 .EQ. 1) GOTO 841 ! yes we have!
    265234      zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
    266       zsat=(local_q(i,k1-1)-v_qs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1))   &
    267         + (local_q(i,k1)-v_qs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
     235      zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1))   &
     236        + (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
    268237      IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
    269238
    … …  
    271240      k1 = k1 - 1
    272241      IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
    273       zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1))               &
     242      zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1))               &
    274243                  *(pplev(i,k1-1)-pplev(i,k1))
    275244      zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
    … …  
    299268      DO i = 1, ngridmx
    300269         IF (itest(i)) THEN
    301          delta_q(i,k) = local_q(i,k) - q(i,k)
     270         delta_q(i,k) = local_q(i,k) - zq(i,k)
    302271         IF (delta_q(i,k).LT.0.) rneb(i,k)  = 1.0
    303272         ENDIF
    … …  
    327296      DO i = 1, ngridmx
    328297         IF (itest(i)) THEN
    329          IF (zq2(i).NE.0.0) &
    330            d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
     298           IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
    331299         ENDIF
    332300      ENDDO
    … …  
    343311      DO k = 1, nlayermx
    344312      DO i = 1, ngridmx
    345          d_t(i,k) = local_t(i,k) - t(i,k)
    346          d_q(i,k) = local_q(i,k) - q(i,k)
     313         d_t(i,k) = local_t(i,k) - zt(i,k)
     314         d_q(i,k) = local_q(i,k) - zq(i,k)
    347315      ENDDO
    348316      ENDDO
    … …  
    352320         DO i = 1, ngridmx
    353321           
    354             dtmana(i,k)       = d_t(i,k)/ptimestep
    355             dqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
    356             dqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
     322            pdtmana(i,k)       = d_t(i,k)/ptimestep
     323            pdqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
     324            pdqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
    357325         
    358326         ENDDO
    359327      ENDDO
    360328
    361 !      print*,'IN MANABE:'
    362 !      print*,'pplev=',pplev
    363 !      print*,'t=',t
    364 !      print*,'d_t=',d_t
    365 !      print*,'d_q=',d_q
    366 !      print*,'local_q=',local_q
    367 !      print*,'q=',q
    368 !      print*,'v_qs(i,k)=',v_qs
    369 !      print*,'v_qsd(i,k)=',v_qsd
    370 !      print*,'cp_delta_t(k)=',cp_delta_t
    371 
    372 !      print*,'d_ql=',d_ql
    373 !      print*,'delta_q=',delta_q
    374 !      print*,'zq1=',zq1
    375 !      print*,'zq2=',zq2
    376 !!      print*,'i_h2o=',i_h2o
    377 !      print*,'i_ice=',i_ice
    378 !
    379 !      print*,'IN MANABE:'
    380 !      print*,'d_q=',d_q
    381 !      print*,'d_ql=',d_ql
    382 !      print*,'dtmana=',d_t
    383 !     stop
    384 !      print*,'gamcpdz at end=',gamcpdz
    385       !  stop   
    386 
    387 !     Some conservation diagnostics...
    388 !      dEtot=0.0
    389 !      dL1tot=0.0
    390 !      dL2tot=0.0
    391 !      dqtot=0.0
    392 !      masse=0.0
    393 !      DO k = 1, nlayermx
    394 !         DO i = 1, ngridmx
    395 !
    396 !            masse = (pplev(i,k) - pplev(i,k+1))/g
    397 !
    398 !            dEtot  = dEtot  + cpp*d_t(i,k)*masse
    399 !            dL1tot = dL1tot + RLVTT*d_ql(i,k)*masse
    400 !            dL2tot = dL2tot + RLVTT*d_q(i,k)*masse ! is this line necessary?
    401 !
    402 !            dqtot = dqtot + (d_q(i,k) + d_ql(i,k))*masse
    403 !
    404 !         ENDDO
    405 !      ENDDO
    406 
    407 !        print*,'In manabe energy change=',dEtot
    408 !        print*,'In manabe condense energy change 1 =',dL1tot
    409 !        print*,'In manabe condense energy change 2 =',dL2tot
    410 !        print*,'In manabe water change=',dqtot
    411329
    412330      RETURN

  • trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F

    r699 r728  
    1111#include "dimensions.h"
    1212#include "dimphys.h"
    13 !#include "comgeomfi.h"
     13#include "comgeomfi.h"
    1414#include "surfdat.h"
    1515#include "planete.h"
    … …  
    111111     .              p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
    112112c
    113 c Read latitudes (coordinates): No need, these are provided by the dynamics
    114 c
    115 !      ierr=nf90_inq_varid(nid,"latitude",nvarid)
    116 !      IF (ierr.NE.nf90_noerr) THEN
    117 !         PRINT*, 'phyetat0: Le champ <latitude> est absent'
    118 !         write(*,*)trim(nf90_strerror(ierr))
    119 !         CALL abort
    120 !      ENDIF
    121 !      ierr=nf90_get_var(nid,nvarid,lati)
    122 !      IF (ierr.NE.nf90_noerr) THEN
    123 !         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
    124 !         write(*,*)trim(nf90_strerror(ierr))
    125 !         CALL abort
    126 !      ENDIF
    127 c
    128 c read longitudes (coordinates): No need, these are provided by the dynamics
    129 c
    130 !      ierr=nf90_inq_varid(nid,"longitude",nvarid)
    131 !      IF (ierr.NE.nf90_noerr) THEN
    132 !         PRINT*, 'phyetat0: Le champ <longitude> est absent'
    133 !         write(*,*)trim(nf90_strerror(ierr))
    134 !         CALL abort
    135 !      ENDIF
    136 !      ierr=nf90_get_var(nid,nvarid,long)
    137 !      IF (ierr.NE.nf90_noerr) THEN
    138 !         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
    139 !         write(*,*)trim(nf90_strerror(ierr))
    140 !         CALL abort
    141 !      ENDIF
    142 c
    143 c Read areas of meshes: No need, these are provided by the dynamics
    144 c
    145 !      ierr=nf90_inq_varid(nid,"area",nvarid)
    146 !      IF (ierr.NE.nf90_noerr) THEN
    147 !         PRINT*, 'phyetat0: Le champ <area> est absent'
    148 !         write(*,*)trim(nf90_strerror(ierr))
    149 !         CALL abort
    150 !      ENDIF
    151 !      ierr=nf90_get_var(nid,nvarid,area)
    152 !      IF (ierr.NE.nf90_noerr) THEN
    153 !         PRINT*, 'phyetat0: Lecture echouee pour <area>'
    154 !         write(*,*)trim(nf90_strerror(ierr))
    155 !         CALL abort
    156 !      ENDIF
    157 !      xmin = 1.0E+20
    158 !      xmax = -1.0E+20
    159 !      xmin = MINVAL(area)
    160 !      xmax = MAXVAL(area)
    161 !      PRINT*,'Aires des mailles <area>:', xmin, xmax
     113c Lecture des latitudes (coordonnees):
     114c
     115      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
     116      IF (ierr.NE.NF_NOERR) THEN
     117         PRINT*, 'phyetat0: Le champ <latitude> est absent'
     118         CALL abort
     119      ENDIF
     120#ifdef NC_DOUBLE
     121      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lati)
     122#else
     123      ierr = NF_GET_VAR_REAL(nid, nvarid, lati)
     124#endif
     125      IF (ierr.NE.NF_NOERR) THEN
     126         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
     127         CALL abort
     128      ENDIF
     129c
     130c Lecture des longitudes (coordonnees):
     131c
     132      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
     133      IF (ierr.NE.NF_NOERR) THEN
     134         PRINT*, 'phyetat0: Le champ <longitude> est absent'
     135         CALL abort
     136      ENDIF
     137#ifdef NC_DOUBLE
     138      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, long)
     139#else
     140      ierr = NF_GET_VAR_REAL(nid, nvarid, long)
     141#endif
     142      IF (ierr.NE.NF_NOERR) THEN
     143         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
     144         CALL abort
     145      ENDIF
     146c
     147c Lecture des aires des mailles:
     148c
     149      ierr = NF_INQ_VARID (nid, "area", nvarid)
     150      IF (ierr.NE.NF_NOERR) THEN
     151         PRINT*, 'phyetat0: Le champ <area> est absent'
     152         CALL abort
     153      ENDIF
     154#ifdef NC_DOUBLE
     155      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, area)
     156#else
     157      ierr = NF_GET_VAR_REAL(nid, nvarid, area)
     158#endif
     159      IF (ierr.NE.NF_NOERR) THEN
     160         PRINT*, 'phyetat0: Lecture echouee pour <area>'
     161         CALL abort
     162      ENDIF
     163      xmin = 1.0E+20
     164      xmax = -1.0E+20
     165      xmin = MINVAL(area)
     166      xmax = MAXVAL(area)
     167      PRINT*,'Aires des mailles <area>:', xmin, xmax
    162168c
    163169c Lecture du geopotentiel au sol:

  • trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

    r726 r728  
    88 
    99      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
    10       use watercommon_h, only : RLVTT
     10      use watercommon_h, only : RLVTT, Psat_water
    1111      use gases_h
    1212      use radcommon_h, only: sigma
     13      use radii_mod, only: h2o_reffrad, co2_reffrad
    1314      use aerosol_mod
    1415      implicit none
    … …  
    191192
    192193      character*80 fichier
    193       integer l,ig,ierr,iq,i, tapphys,nw
     194      integer l,ig,ierr,iq,iaer
    194195
    195196      real fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
    … …  
    213214      real zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
    214215      real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
    215 
    216       real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
    217216
    218217!     Tendencies due to various processes:
    … …  
    231230      real zdhadj(ngridmx,nlayermx)                           ! (K/s)
    232231      real zdtgw(ngridmx,nlayermx)                            ! (K/s)
     232      real zdtmr(ngridmx,nlayermx)
    233233      real zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
    234234      real zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
    235235      real zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
     236      real zdumr(ngridmx,nlayermx),zdvmr(ngridmx,nlayermx)
     237      real zdtsurfmr(ngridmx)
     238     
     239      real zdmassmr(ngridmx,nlayermx),zdpsrfmr(ngridmx)
    236240
    237241      real zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
     242      real zdqevap(ngridmx,nlayermx)
    238243      real zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
    239244      real zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
    240245      real zdqadj(ngridmx,nlayermx,nqmx)
    241246      real zdqc(ngridmx,nlayermx,nqmx)
     247      real zdqmr(ngridmx,nlayermx,nqmx),zdqsurfmr(ngridmx,nqmx)
    242248      real zdqlscale(ngridmx,nlayermx,nqmx)
    243249      real zdqslscale(ngridmx,nqmx)
    … …  
    329335      real qevap(ngridmx,nlayermx,nqmx)
    330336      real tevap(ngridmx,nlayermx)
    331       real dqevap(ngridmx,nlayermx)
    332       real dtevap(ngridmx,nlayermx)
     337      real dqevap1(ngridmx,nlayermx)
     338      real dtevap1(ngridmx,nlayermx)
    333339
    334340!     included by BC for hydrology
    … …  
    343349
    344350!     included by RW for RH diagnostic
    345       real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx)
     351      real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx),psat_tmp
    346352
    347353!     included by RW for hydrology
    … …  
    366372
    367373!     included by BC for cloud fraction computations
    368       real cloudfrac(ngridmx,nlayermx)
    369       real totcloudfrac(ngridmx)
     374      real,save :: cloudfrac(ngridmx,nlayermx)
     375      real,save :: totcloudfrac(ngridmx)
    370376
    371377!     included by RW for vdifc water conservation test
    … …  
    389395
    390396!     included by RW for variable H2O particle sizes
     397      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
     398      real :: reffrad1(ngridmx,nlayermx,naerkind) ! for clearsky call to callcorrk
     399      real,save :: nueffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius variance
     400      real :: reffh2oliq(ngridmx,nlayermx) ! liquid water particles effective radius (m)
     401      real :: reffh2oice(ngridmx,nlayermx) ! water ice particles effective radius (m)
    391402      real reffH2O(ngridmx,nlayermx)
    392403      real reffcol(ngridmx,naerkind)
    393404
    394405!     included by RW for sourceevol
    395       real ice_initial(ngridmx)!, isoil
     406      real, save :: ice_initial(ngridmx)!, isoil
    396407      real delta_ice,ice_tot
    397       real ice_min(ngridmx)
    398       save ice_min
    399       save ice_initial
    400 
     408      real, save :: ice_min(ngridmx)
     409     
    401410      integer num_run
    402       logical ice_update
    403       save ice_update
    404 
     411      logical,save :: ice_update
     412     
    405413!=======================================================================
    406414     
    … …  
    457465         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
    458466
    459          do ig=1,ngrid
    460             albedo(ig)=albedo0(ig)
    461          enddo
     467         albedo(:)=albedo0(:)
    462468
    463469         if(tlocked)then
    … …  
    473479         else
    474480            print*,'WARNING! Thermal conduction in the soil turned off'
    475             do ig=1,ngrid
    476                capcal(ig)=1.e6
    477                fluxgrd(ig)=0.
    478                if(noradsurf)then
    479                   fluxgrd(ig)=10.0
    480                endif
    481             enddo
     481            capcal(:)=1.e6
     482            fluxgrd(:)=0.
     483            if(noradsurf)then
     484                  fluxgrd(:)=10.0
     485            endif
    482486            print*,'Flux from ground = ',fluxgrd,' W m^-2'
    483487         endif
    … …  
    502506!        define surface as continent or ocean
    503507!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     508         rnat(:)=1
    504509         do ig=1,ngridmx
    505             rnat(ig)=1
    506 
    507510!            if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
    508511            if(inertiedat(ig,1).gt.1.E4)then
    … …  
    517520!        initialise surface history variable
    518521!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    519          do ig=1,ngridmx
    520             do iq=1,nqmx
    521                qsurf_hist(ig,iq)=qsurf(ig,iq)
    522             enddo
    523          enddo
     522         qsurf_hist(:,:)=qsurf(:,:)
    524523
    525524!        initialise variable for dynamical heating diagnostic
    … …  
    556555            endif
    557556
    558             do ig=1,ngrid
    559 !already done above (JL12)
    560 !               qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
    561                if(ice_update)then
    562                   ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
    563                endif
    564             enddo
     557            if(ice_update)then
     558               ice_initial(:)=qsurf(:,igcm_h2o_ice)
     559            endif
    565560
    566561         endif
    … …  
    585580!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    586581
    587       pdu(:,:) = 0.0
    588       pdv(:,:) = 0.0
    589 !      if ( (.not.nearco2cond).and.(.not.firstcall) ) then
     582      pdu(1:ngridmx,1:nlayermx) = 0.0
     583      pdv(1:ngridmx,1:nlayermx) = 0.0
    590584      if ( .not.nearco2cond ) then
    591          pdt(:,:) = 0.0
    592       endif ! this was the source of an evil bug...
    593       pdq(:,:,:) = 0.0
    594       pdpsrf(:)   = 0.0
    595       zflubid(:)  = 0.0
    596       zdtsurf(:)  = 0.0
    597       dqsurf(:,:) = 0.0
     585         pdt(1:ngridmx,1:nlayermx) = 0.0
     586      endif
     587      pdq(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
     588      pdpsrf(1:ngridmx)       = 0.0
     589      zflubid(1:ngridmx)      = 0.0
     590      zdtsurf(1:ngridmx)      = 0.0
     591      dqsurf(1:ngridmx,1:nqmx)= 0.0
    598592
    599593      zday=pday+ptime ! compute time, in sols (and fraction thereof)
    … …  
    610604!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    611605
    612       do l=1,nlayer
    613          do ig=1,ngrid
    614             zzlay(ig,l)=pphi(ig,l)/g
    615          enddo
    616       enddo
    617       do ig=1,ngrid
    618          zzlev(ig,1)=0.
    619          zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
    620       enddo
     606      zzlay(1:ngridmx,1:nlayermx)=pphi(1:ngridmx,1:nlayermx)/g
     607
     608      zzlev(1:ngridmx,1)=0.
     609      zzlev(1:ngridmx,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
     610
    621611      do l=2,nlayer
    622612         do ig=1,ngrid
    … …  
    631621!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    632622
    633 
    634623      do l=1,nlayer         
    635624         do ig=1,ngrid
    … …  
    657646              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
    658647
    659                call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
     648              call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
    660649               ztim1,ztim2,ztim3,mu0,fract)
    661650
    … …  
    707696                      zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
    708697                      fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
    709                       reffrad,tau_col1,cloudfrac,totcloudfrac,             &
     698                      reffrad1,tau_col1,cloudfrac,totcloudfrac,             &
    710699                      clearsky,firstcall,lastcall)
    711700                 clearsky = .false.  ! just in case.     
    … …  
    722711                    tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
    723712                   
    724                     do l=1,nlayer
    725                        zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw(ig,l)
    726                        zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw(ig,l)
    727                     enddo
    728 
    729                     do nw=1,L_NSPECTV
    730                        OSR_nu(ig,nw) = ntf*OSR_nu1(ig,nw) + tf*OSR_nu(ig,nw)                   
    731                     enddo
    732                     do nw=1,L_NSPECTI
    733                        OLR_nu(ig,nw) = ntf*OLR_nu1(ig,nw) + tf*OLR_nu(ig,nw)                   
    734                     enddo
     713                    zdtlw(ig,1:nlayermx) = ntf*zdtlw1(ig,1:nlayermx) + tf*zdtlw(ig,1:nlayermx)
     714                    zdtsw(ig,1:nlayermx) = ntf*zdtsw1(ig,1:nlayermx) + tf*zdtsw(ig,1:nlayermx)
     715
     716                    OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)                   
     717                    OLR_nu(ig,1:L_NSPECTV) = ntf*OLR_nu1(ig,1:L_NSPECTV) + tf*OLR_nu(ig,1:L_NSPECTV)                   
    735718
    736719                 enddo
    … …  
    740723!             Radiative flux from the sky absorbed by the surface (W.m-2)
    741724              GSR=0.0
    742               do ig=1,ngrid
    743                  fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)             &
    744                       +fluxsurf_sw(ig)*(1.-albedo(ig))
    745 
    746                  if(noradsurf)then ! no lower surface; SW flux just disappears
    747                     GSR = GSR + fluxsurf_sw(ig)*area(ig)
    748                     fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
    749                  endif
    750 
    751               enddo
    752               if(noradsurf)then
    753                  print*,'SW lost in deep atmosphere = ',GSR/totarea,' W m^-2'
     725              fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)+fluxsurf_sw(1:ngridmx)*(1.-albedo(1:ngridmx))
     726
     727              if(noradsurf)then ! no lower surface; SW flux just disappears
     728                 GSR = SUM(fluxsurf_sw(1:ngridmx)*area(1:ngridmx))/totarea
     729                 fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)
     730                 print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
    754731              endif
    755732
    756733!             Net atmospheric radiative heating rate (K.s-1)
    757               do l=1,nlayer
    758                  do ig=1,ngrid
    759                     dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
    760                  enddo
    761               enddo
    762 
     734              dtrad(1:ngridmx,1:nlayermx)=zdtsw(1:ngridmx,1:nlayermx)+zdtlw(1:ngridmx,1:nlayermx)
    763735
    764736           elseif(newtonian)then
    … …  
    768740              call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
    769741
    770               do ig=1,ngrid
    771                  zdtsurf(ig) = +(pt(ig,1)-tsurf(ig))/ptimestep
    772               enddo
     742              zdtsurf(1:ngridmx) = +(pt(1:ngridmx,1)-tsurf(1:ngridmx))/ptimestep
    773743              ! e.g. surface becomes proxy for 1st atmospheric layer ?
    774744
    … …  
    777747!          c) Atmosphere has no radiative effect
    778748!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    779               do ig=1,ngrid
    780                  fluxtop_dn(ig)  = fract(ig)*mu0(ig)*Fat1AU/dist_star**2
    781                  if(ngrid.eq.1)then ! / by 4 globally in 1D case...
    782                     fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
    783                  endif
    784                  fluxsurf_sw(ig) = fluxtop_dn(ig)
    785                  fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
    786                  fluxtop_lw(ig)  = emis(ig)*sigma*tsurf(ig)**4
    787               enddo ! radiation skips the atmosphere entirely
    788 
    789               do l=1,nlayer
    790                  do ig=1,ngrid
    791                     dtrad(ig,l)=0.0
    792                  enddo
    793               enddo ! hence no atmospheric radiative heating
     749              fluxtop_dn(1:ngridmx)  = fract(1:ngridmx)*mu0(1:ngridmx)*Fat1AU/dist_star**2
     750              if(ngrid.eq.1)then ! / by 4 globally in 1D case...
     751                 fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
     752              endif
     753              fluxsurf_sw(1:ngridmx) = fluxtop_dn(1:ngridmx)
     754              fluxrad_sky(1:ngridmx) = fluxtop_dn(1:ngridmx)*(1.-albedo(1:ngridmx))
     755              fluxtop_lw(1:ngridmx)  = emis(1:ngridmx)*sigma*tsurf(1:ngridmx)**4
     756              ! radiation skips the atmosphere entirely
     757
     758
     759              dtrad(1:ngridmx,1:nlayermx)=0.0
     760              ! hence no atmospheric radiative heating
    794761
    795762           endif                ! if corrk
    … …  
    801768!    ------------------------------------------
    802769
    803         do ig=1,ngrid
    804            zplanck(ig)=tsurf(ig)*tsurf(ig)
    805            zplanck(ig)=emis(ig)*sigma*zplanck(ig)*zplanck(ig)
    806            fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
    807         enddo
    808 
    809         do l=1,nlayer
    810            do ig=1,ngrid
    811               pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
    812            enddo
    813         enddo
     770        zplanck(1:ngridmx)=tsurf(1:ngridmx)*tsurf(1:ngridmx)
     771        zplanck(1:ngridmx)=emis(1:ngridmx)*sigma*zplanck(1:ngridmx)*zplanck(1:ngridmx)
     772        fluxrad(1:ngridmx)=fluxrad_sky(1:ngridmx)-zplanck(1:ngridmx)
     773        pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+dtrad(1:ngridmx,1:nlayermx)
    814774
    815775!-------------------------
    … …  
    840800      if (calldifv) then
    841801     
    842          do ig=1,ngrid
    843             zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
    844          enddo
    845 
    846          zdum1(:,:)=0.0
    847          zdum2(:,:)=0.0
     802         zflubid(1:ngridmx)=fluxrad(1:ngridmx)+fluxgrd(1:ngridmx)
     803
     804         zdum1(1:ngridmx,1:nlayermx)=0.0
     805         zdum2(1:ngridmx,1:nlayermx)=0.0
    848806
    849807
    … …  
    857815              zdum1,zdum2,pdt,pdq,zflubid,           &
    858816              zdudif,zdvdif,zdtdif,zdtsdif,          &
    859               sensibFlux,q2,zdqdif,zdqsdif,lastcall)
     817              sensibFlux,q2,zdqdif,zdqevap,zdqsdif,lastcall)
    860818
    861819         else
    862820         
    863            do l=1,nlayer
    864              do ig=1,ngrid
    865                zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
    866              enddo
    867            enddo
     821           zdh(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
    868822 
    869823           call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
    … …  
    875829              sensibFlux,q2,zdqdif,zdqsdif,lastcall)
    876830
    877            do l=1,nlayer
    878              do ig=1,ngrid
    879                 zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
    880              enddo
    881            enddo
     831           zdtdif(1:ngridmx,1:nlayermx)=zdhdif(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
     832           zdqevap(1:ngridmx,1:nlayermx)=0.
    882833
    883834         end if !turbdiff
    884835
    885          do l=1,nlayer
    886             do ig=1,ngrid
    887                pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
    888                pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
    889                pdt(ig,l)=pdt(ig,l)+zdtdif(ig,l)
    890             enddo
    891          enddo
    892 
    893          do ig=1,ngrid
    894             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
    895          enddo
    896 
     836         pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvdif(1:ngridmx,1:nlayermx)
     837         pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zdudif(1:ngridmx,1:nlayermx)
     838         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtdif(1:ngridmx,1:nlayermx)
     839         zdtsurf(1:ngridmx)=zdtsurf(1:ngridmx)+zdtsdif(1:ngridmx)
    897840         if (tracer) then
    898            do iq=1, nq
    899             do l=1,nlayer
    900               do ig=1,ngrid
    901                  pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
    902               enddo
    903             enddo
    904            enddo
    905            do iq=1, nq
    906               do ig=1,ngrid
    907                  dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
    908               enddo
    909            enddo
    910 
     841           pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqdif(1:ngridmx,1:nlayermx,1:nqmx)
     842           dqsurf(1:ngridmx,1:nqmx)=dqsurf(1:ngridmx,1:nqmx) + zdqsdif(1:ngridmx,1:nqmx)
    911843         end if ! of if (tracer)
    912844
    … …  
    965897         if(.not.newtonian)then
    966898
    967             do ig=1,ngrid
    968                zdtsurf(ig) = zdtsurf(ig) + (fluxrad(ig) + fluxgrd(ig))/capcal(ig)
    969             enddo
     899            zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + (fluxrad(1:ngridmx) + fluxgrd(1:ngridmx))/capcal(1:ngridmx)
    970900
    971901         endif
    … …  
    980910      if(calladj) then
    981911
    982          do l=1,nlayer
    983             do ig=1,ngrid
    984                zdh(ig,l) = pdt(ig,l)/zpopsk(ig,l)
    985             enddo
    986          enddo
    987          zduadj(:,:)=0.0
    988          zdvadj(:,:)=0.0
    989          zdhadj(:,:)=0.0
     912         zdh(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
     913         zduadj(1:ngridmx,1:nlayermx)=0.0
     914         zdvadj(1:ngridmx,1:nlayermx)=0.0
     915         zdhadj(1:ngridmx,1:nlayermx)=0.0
    990916
    991917
    … …  
    997923              zdqadj)
    998924
    999          do l=1,nlayer
    1000             do ig=1,ngrid
    1001                pdu(ig,l) = pdu(ig,l) + zduadj(ig,l)
    1002                pdv(ig,l) = pdv(ig,l) + zdvadj(ig,l)
    1003                pdt(ig,l)    = pdt(ig,l) + zdhadj(ig,l)*zpopsk(ig,l)
    1004                zdtadj(ig,l) = zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
    1005             enddo
    1006          enddo
    1007 
     925         pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zduadj(1:ngridmx,1:nlayermx)
     926         pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvadj(1:ngridmx,1:nlayermx)
     927         pdt(1:ngridmx,1:nlayermx)    = pdt(1:ngridmx,1:nlayermx) + zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx)
     928         zdtadj(1:ngridmx,1:nlayermx) = zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
     929 
    1008930         if(tracer) then
    1009            do iq=1, nq
    1010             do l=1,nlayer
    1011               do ig=1,ngrid
    1012                  pdq(ig,l,iq) = pdq(ig,l,iq) + zdqadj(ig,l,iq)
    1013               enddo
    1014             enddo
    1015            enddo
     931            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqadj(1:ngridmx,1:nlayermx,1:nqmx)
    1016932         end if
    1017933
    … …  
    1061977              zdqc,reffrad)
    1062978
    1063          do l=1,nlayer
    1064            do ig=1,ngrid
    1065              pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
    1066              pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
    1067              pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
    1068            enddo
    1069          enddo
    1070          do ig=1,ngrid
    1071             zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
    1072          enddo
    1073 
    1074          do iq=1,nq ! should use new notation here !
    1075             do l=1,nlayer
    1076                do ig=1,ngrid
    1077                   pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
    1078                enddo
    1079             enddo
    1080          enddo
     979
     980         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtc(1:ngridmx,1:nlayermx)
     981         pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvc(1:ngridmx,1:nlayermx)
     982         pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zduc(1:ngridmx,1:nlayermx)
     983         zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfc(1:ngridmx)
     984
     985         pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqc(1:ngridmx,1:nlayermx,1:nqmx)
    1081986         ! Note: we do not add surface co2ice tendency
    1082987         ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
    … …  
    11081013!        Water ice condensation in the atmosphere
    11091014!        ----------------------------------------
    1110             if(watercond)then
    1111 
    1112                if(RLVTT.gt.1.e-8)then
    1113 
     1015            if(watercond.and.(RLVTT.gt.1.e-8))then
     1016
     1017!             ----------------
     1018!             Moist convection
     1019!             ----------------
    11141020               ! Re-evaporate cloud water/ice
    1115                call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
    1116                DO l = 1, nlayer 
    1117                   DO ig = 1, ngrid
    1118                      pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
    1119                      pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
    1120                      pdt(ig,l)              = pdt(ig,l)+dtevap(ig,l)
    1121                   enddo
    1122                enddo
    1123 
    1124                call moistadj(pt,qevap,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
    1125                do l=1,nlayer
    1126                   do ig=1,ngrid
    1127                      pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqmoist(ig,l,igcm_h2o_vap)
    1128                      pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqmoist(ig,l,igcm_h2o_ice)
    1129                      pdt(ig,l)              = pdt(ig,l)+dtmoist(ig,l)
    1130                   enddo
    1131                enddo
     1021!              dqevap1(1:ngridmx,1:nlayermx)=0.
     1022!              dtevap1(1:ngridmx,1:nlayermx)=0.
     1023!               call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap)
     1024               
     1025               dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
     1026               dtmoist(1:ngridmx,1:nlayermx)=0.
     1027
     1028               call moistadj(pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
     1029
     1030               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)   &
     1031                          +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)
     1032               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) =pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)     &
     1033                          +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
     1034               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
     1035
     1036!               do l=1,nlayermx
     1037!                do ig=1,ngridmx
     1038!                 if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l
     1039!                 if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l
     1040!                 if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l               
     1041!                  if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP
     1042!                enddo
     1043!              enddo
     1044!              print*,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)*86400.
    11321045
    11331046               !-------------------------
    11341047               ! test energy conservation
    11351048               if(enertest)then
    1136                   dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
     1049                  dEtot=cpp*SUM(massarea(:,:)*dtmoist(:,:))/totarea
    11371050                  do ig = 1, ngrid
    1138                      madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
     1051                     madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
    11391052                  enddo
    11401053                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
     1054                  print*,'In moistadj MAX atmospheric energy change   =',MAXVAL(pdt(:,:))*86400.,'K/day'
    11411055                 
    11421056                ! test energy conservation
    … …  
    11461060               endif
    11471061               !-------------------------
    1148                
    1149 
    1150                endif
    11511062               
    11521063
    1153                ! Re-evaporate cloud water/ice
    1154                call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
    1155                do l = 1, nlayer 
    1156                   do ig = 1, ngrid
    1157                      pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
    1158                      pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
    1159                      pdt(ig,l)              = pdt(ig,l)+dtevap(ig,l)
    1160                   enddo
    1161                enddo ! note: we use qevap but not tevap in largescale/moistadj
    1162                      ! otherwise is a big mess
    1163 
    1164                call largescale(ptimestep,pplev,pplay,pt,qevap,        & ! a bug was here!
    1165                     pdt,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc,reffH2O)
    1166                do l=1,nlayer
    1167                   do ig=1,ngrid
    1168                      pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqvaplscale(ig,l)
    1169                      pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l)
    1170                      pdt(ig,l)              = pdt(ig,l)+dtlscale(ig,l)
    1171                      
    1172                      if(aeroh2o.and.(.not.aerofixh2o)) then
    1173                         reffrad(ig,l,iaero_h2o)=reffH2O(ig,l)
    1174                      endif
    1175 
    1176                   enddo
    1177                enddo
     1064!        --------------------------------
     1065!        Large scale condensation/evaporation
     1066!        --------------------------------
     1067
     1068               call largescale(ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
     1069
     1070               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtlscale(1:ngridmx,1:nlayermx)
     1071               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngridmx,1:nlayermx)
     1072               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngridmx,1:nlayermx)
    11781073
    11791074               !-------------------------
    11801075               ! test energy conservation
    11811076               if(enertest)then
    1182                   dEtot=cpp*SUM(massarea(:,:)*(dtmoist(:,:)+dtevap(:,:)))/totarea
    11831077                  do ig = 1, ngrid
    1184                      madjdE(ig) = cpp*SUM(mass(:,:)*(dtmoist(:,:)+dtevap(:,:)))
     1078                     lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
    11851079                  enddo
    1186                   print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
    1187                  
    1188                 ! test energy conservation
    1189                   dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
    1190                   dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
    1191                   print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
     1080                  dEtot=cpp*SUM(massarea(:,:)*(dtlscale(:,:)))/totarea
     1081                  if(isnan(dEtot)) then
     1082                     print*,'Nan in largescale, abort'
     1083                     STOP
     1084                  endif
     1085                  print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
     1086
     1087               ! test water conservation
     1088                  dWtot = SUM(massarea(:,:)*dqvaplscale(:,:))*ptimestep/totarea
     1089                  dWtot = dWtot + SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea
     1090                  print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
    11921091               endif
    11931092               !-------------------------
    … …  
    12001099               enddo
    12011100
    1202                ! compute total cloud fraction in column
    1203                call totalcloudfrac(cloudfrac,totcloudfrac)
    1204 
    1205            endif                ! of if (watercondense)
     1101            endif                ! of if (watercondense)
    12061102           
    12071103
    … …  
    12091105!        Water ice / liquid precipitation
    12101106!        --------------------------------
    1211            if(waterrain)then
    1212 
    1213               zdqrain(:,:,:) = 0.0
    1214               zdqsrain(:)    = 0.0
    1215               zdqssnow(:)    = 0.0
    1216 
    1217               call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
     1107            if(waterrain)then
     1108
     1109               zdqrain(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
     1110               zdqsrain(1:ngridmx)    = 0.0
     1111               zdqssnow(1:ngridmx)    = 0.0
     1112
     1113               call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
    12181114                   zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
    12191115
    1220               do l=1,nlayer
    1221                  do ig=1,ngrid
    1222                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+zdqrain(ig,l,igcm_h2o_vap)
    1223                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+zdqrain(ig,l,igcm_h2o_ice)
    1224                     pdt(ig,l)              = pdt(ig,l)+zdtrain(ig,l)
    1225                  enddo
    1226               enddo
    1227 
    1228               do ig=1,ngrid
    1229                  dqsurf(ig,igcm_h2o_vap) = dqsurf(ig,igcm_h2o_vap)+zdqsrain(ig) ! a bug was here
    1230                  dqsurf(ig,igcm_h2o_ice) = dqsurf(ig,igcm_h2o_ice)+zdqssnow(ig)
    1231                  rainout(ig)             = zdqsrain(ig)+zdqssnow(ig) ! diagnostic   
    1232               enddo
    1233 
     1116               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
     1117                     +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)
     1118               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
     1119                     +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_ice)
     1120               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+zdtrain(1:ngridmx,1:nlayermx)
     1121               dqsurf(1:ngridmx,igcm_h2o_vap) = dqsurf(1:ngridmx,igcm_h2o_vap)+zdqsrain(1:ngridmx) ! a bug was here
     1122               dqsurf(1:ngridmx,igcm_h2o_ice) = dqsurf(1:ngridmx,igcm_h2o_ice)+zdqssnow(1:ngridmx)
     1123               rainout(1:ngridmx)             = zdqsrain(1:ngridmx)+zdqssnow(1:ngridmx) ! diagnostic   
    12341124
    12351125
    … …  
    12421132                  dItot = dItot + SUM(area(:)*zdqssnow(:))/totarea*RLVTT/cpp
    12431133                  dVtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
    1244                   dVtot = dItot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
     1134                  dVtot = dVtot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
    12451135                  dEtot = dItot + dVtot
    12461136                 print*,'In rain dItot =',dItot,' W m-2'
    … …  
    12581148              !-------------------------
    12591149
    1260            end if                 ! of if (waterrain)
    1261         end if                    ! of if (water)
     1150            end if                 ! of if (waterrain)
     1151         end if                    ! of if (water)
    12621152
    12631153
    … …  
    12681158!        -------------
    12691159        if (sedimentation) then
    1270            zdqsed(:,:,:) = 0.0
    1271            zdqssed(:,:)  = 0.0
     1160           zdqsed(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
     1161           zdqssed(1:ngridmx,1:nqmx)  = 0.0
    12721162
    12731163
    … …  
    12971187           !-------------------------
    12981188
    1299            do iq=1,nq
    1300               ! for now, we only allow H2O ice to sediment
     1189           ! for now, we only allow H2O ice to sediment
    13011190              ! and as in rain.F90, whether it falls as rain or snow depends
    13021191              ! only on the surface temperature
    1303               do ig=1,ngrid
    1304                  do l=1,nlayer
    1305                     pdq(ig,l,iq) = pdq(ig,l,iq) + zdqsed(ig,l,iq)
    1306                  enddo
    1307                  dqsurf(ig,iq) = dqsurf(ig,iq) + zdqssed(ig,iq)
    1308               enddo
    1309            enddo
     1192           pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqsed(1:ngridmx,1:nlayermx,1:nqmx)
     1193           dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqssed(1:ngridmx,1:nqmx)
    13101194
    13111195           !-------------------------
    … …  
    13261210!     ---------
    13271211
     1212!       -----------------------------------
     1213!       Updating atm mass and tracer budget
     1214!       -----------------------------------
     1215
     1216         if(mass_redistrib) then
     1217
     1218            zdmassmr(1:ngridmx,1:nlayermx) = mass(1:ngridmx,1:nlayermx) * &
     1219               (   zdqevap(1:ngridmx,1:nlayermx)                          &
     1220                 + zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
     1221                 + dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
     1222                 + dqvaplscale(1:ngridmx,1:nlayermx) )
     1223                 
     1224           
     1225            call writediagfi(ngridmx,"mass_evap","mass gain"," ",3,zdmassmr)
     1226            call writediagfi(ngridmx,"mass","mass"," ",3,mass)
     1227
     1228            call mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
     1229                       rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
     1230                       pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,  &
     1231                       zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
     1232         
     1233
     1234            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqmr(1:ngridmx,1:nlayermx,1:nqmx)
     1235            dqsurf(1:ngridmx,1:nqmx)         = dqsurf(1:ngridmx,1:nqmx) + zdqsurfmr(1:ngridmx,1:nqmx)
     1236            pdt(1:ngridmx,1:nlayermx)        = pdt(1:ngridmx,1:nlayermx) + zdtmr(1:ngridmx,1:nlayermx)
     1237            pdu(1:ngridmx,1:nlayermx)        = pdu(1:ngridmx,1:nlayermx) + zdumr(1:ngridmx,1:nlayermx)
     1238            pdv(1:ngridmx,1:nlayermx)        = pdv(1:ngridmx,1:nlayermx) + zdvmr(1:ngridmx,1:nlayermx)
     1239            pdpsrf(1:ngridmx)                = pdpsrf(1:ngridmx) + zdpsrfmr(1:ngridmx)
     1240            zdtsurf(1:ngridmx)               = zdtsurf(1:ngridmx) + zdtsurfmr(1:ngridmx)
     1241           
     1242!           print*,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayermx,igcm_h2o_vap)
     1243         endif
     1244
     1245
     1246
    13281247!       ---------------------------------
    13291248!       Updating tracer budget on surface
    … …  
    13361255            ! note: for now, also changes albedo in the subroutine
    13371256
    1338             do ig=1,ngrid
    1339                zdtsurf(ig) = zdtsurf(ig) + zdtsurf_hyd(ig)
    1340                do iq=1,nq
    1341                   qsurf(ig,iq) = qsurf(ig,iq)+ptimestep*dqs_hyd(ig,iq)
    1342                enddo
    1343             enddo
     1257            zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurf_hyd(1:ngridmx)
     1258            qsurf(1:ngridmx,1:nqmx) = qsurf(1:ngridmx,1:nqmx) + ptimestep*dqs_hyd(1:ngridmx,1:nqmx)
    13441259            ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
    13451260
    … …  
    13651280         ELSE     ! of if (hydrology)
    13661281
    1367             do iq=1,nq
    1368                do ig=1,ngrid
    1369                   qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
    1370                enddo
    1371             enddo           
     1282            qsurf(1:ngridmx,1:nqmx)=qsurf(1:ngridmx,1:nqmx)+ptimestep*dqsurf(1:ngridmx,1:nqmx)
    13721283
    13731284         END IF   ! of if (hydrology)
    … …  
    13801291
    13811292         if(ice_update)then
    1382             do ig = 1, ngrid
    1383                ice_min(ig)=min(ice_min(ig),qsurf(ig,igcm_h2o_ice))
    1384             enddo
     1293            ice_min(1:ngridmx)=min(ice_min(1:ngridmx),qsurf(1:ngridmx,igcm_h2o_ice))
    13851294         endif
    13861295
    … …  
    13931302
    13941303!     Increment surface temperature
    1395       do ig=1,ngrid
    1396          tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
    1397       enddo
     1304      tsurf(1:ngridmx)=tsurf(1:ngridmx)+ptimestep*zdtsurf(1:ngridmx)
    13981305
    13991306!     Compute soil temperatures and subsurface heat flux
    … …  
    14061313! test energy conservation
    14071314      if(enertest)then
    1408          dEtots=0.0
    1409          do ig = 1, ngrid
    1410             dEtots = dEtots + capcal(ig)*zdtsurf(ig)*area(ig)
    1411          enddo
    1412          dEtots=dEtots/totarea
     1315         dEtots = SUM(area(:)*capcal(:)*zdtsurf(:))/totarea
    14131316         print*,'Surface energy change                 =',dEtots,' W m-2'
    14141317      endif
    … …  
    14251328 
    14261329      ! temperature, zonal and meridional wind
    1427       do l=1,nlayer
    1428         do ig=1,ngrid
    1429           zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
    1430           zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
    1431           zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
    1432 
    1433           ! diagnostic
    1434           zdtdyn(ig,l)     = ztprevious(ig,l)-pt(ig,l)
    1435           ztprevious(ig,l) = zt(ig,l)
    1436         enddo
     1330      zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx) + pdt(1:ngridmx,1:nlayermx)*ptimestep
     1331      zu(1:ngridmx,1:nlayermx) = pu(1:ngridmx,1:nlayermx) + pdu(1:ngridmx,1:nlayermx)*ptimestep
     1332      zv(1:ngridmx,1:nlayermx) = pv(1:ngridmx,1:nlayermx) + pdv(1:ngridmx,1:nlayermx)*ptimestep
     1333
     1334      ! diagnostic
     1335      zdtdyn(1:ngridmx,1:nlayermx)     = pt(1:ngridmx,1:nlayermx)-ztprevious(1:ngridmx,1:nlayermx)
     1336      ztprevious(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
     1337
     1338      if(firstcall)then
     1339         zdtdyn(1:ngridmx,1:nlayermx)=0.0
     1340      endif
     1341
     1342      ! dynamical heating diagnostic
     1343      do ig=1,ngrid
     1344         fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp/ptimestep
    14371345      enddo
    14381346
    1439       if(firstcall)then
    1440          zdtdyn(:,:)=0.0
    1441       endif
    1442 
    1443       ! dynamical heating diagnostic
    1444       fluxdyn(:)=0.
    1445       do ig=1,ngrid
    1446          do l=1,nlayer
    1447             fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)   &
    1448                  *(pplev(ig,l)-pplev(ig,l+1))*cpp/g
    1449          enddo
    1450       enddo
    1451 
    14521347      ! tracers
    1453       do iq=1, nq
    1454         do l=1,nlayer
    1455           do ig=1,ngrid
    1456             zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
    1457           enddo
    1458         enddo
    1459       enddo
     1348      zq(1:ngridmx,1:nlayermx,1:nqmx) = pq(1:ngridmx,1:nlayermx,1:nqmx) + pdq(1:ngridmx,1:nlayermx,1:nqmx)*ptimestep
    14601349
    14611350      ! surface pressure
    1462       do ig=1,ngrid
    1463           ps(ig) = pplev(ig,1) + pdpsrf(ig)*ptimestep
    1464       enddo
     1351      ps(1:ngridmx) = pplev(1:ngridmx,1) + pdpsrf(1:ngridmx)*ptimestep
    14651352
    14661353      ! pressure
    14671354      do l=1,nlayer
    1468         do ig=1,ngrid
    1469              zplev(ig,l) = pplev(ig,l)/pplev(ig,1)*ps(ig)
    1470              zplay(ig,l) = pplay(ig,l)/pplev(ig,1)*ps(ig)
    1471         enddo
     1355          zplev(1:ngridmx,l) = pplev(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(:)
     1356          zplay(1:ngridmx,l) = pplay(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(1:ngridmx)
    14721357      enddo
    14731358
    … …  
    15581443!     ---------------------------------------------------------
    15591444      if(tracer)then
    1560          qcol(:,:)=0.0
     1445         qcol(1:ngridmx,1:nqmx)=0.0
    15611446         do iq=1,nq
    1562             do ig=1,ngrid
    1563                do l=1,nlayer
    1564                   qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *    &
    1565                             (pplev(ig,l) - pplev(ig,l+1)) / g
    1566                enddo
     1447           do ig=1,ngridmx
     1448              qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
     1449           enddo
     1450         enddo
     1451
     1452         ! Generalised for arbitrary aerosols now. (LK)
     1453         reffcol(1:ngridmx,1:naerkind)=0.0
     1454         if(co2cond.and.(iaero_co2.ne.0))then
     1455            call co2_reffrad(zq,reffrad)
     1456            do ig=1,ngridmx
     1457               reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)*reffrad(ig,1:nlayermx,iaero_co2)*mass(ig,1:nlayermx))
    15671458            enddo
    1568          enddo
    1569 
    1570          ! Generalised for arbitrary aerosols now. (LK)
    1571          reffcol(:,:)=0.0
    1572          do ig=1,ngrid
    1573             do l=1,nlayer
    1574                if(co2cond.and.(iaero_co2.ne.0)) then
    1575                   reffcol(ig,iaero_co2) = reffcol(ig,iaero_co2) + zq(ig,l,igcm_co2_ice) * &
    1576                                   reffrad(ig,l,iaero_co2) *    &
    1577                                   (pplev(ig,l) - pplev(ig,l+1)) / g
    1578                endif
    1579                if(water.and.(iaero_h2o.ne.0)) then
    1580                   reffcol(ig,iaero_h2o) = reffcol(ig,iaero_h2o) + zq(ig,l,igcm_h2o_ice) * &
    1581                                   reffrad(ig,l,iaero_h2o) *    &
    1582                                   (pplev(ig,l) - pplev(ig,l+1)) / g
    1583                endif
     1459         endif
     1460         if(water.and.(iaero_h2o.ne.0))then
     1461            call h2o_reffrad(zq,zt,reffrad,nueffrad)
     1462            do ig=1,ngridmx
     1463               reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)*reffrad(ig,1:nlayermx,iaero_h2o)*mass(ig,1:nlayermx))
    15841464            enddo
    1585          enddo
     1465         endif
    15861466
    15871467      endif
    … …  
    16221502         do l = 1, nlayer
    16231503            do ig = 1, ngrid
    1624                call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
     1504!               call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
     1505               call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
     1506
    16251507               RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
    16261508            enddo
    … …  
    16291511         ! compute maximum possible H2O column amount (100% saturation)
    16301512         do ig=1,ngrid
    1631             H2Omaxcol(ig)=0.0
    1632             do l=1,nlayer
    1633                H2Omaxcol(ig) = H2Omaxcol(ig) + qsat(ig,l) *    &
    1634                     (pplev(ig,l) - pplev(ig,l+1))/g
    1635             enddo
     1513               H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
    16361514         enddo
    16371515
    … …  
    16621540
    16631541                  ! add multiple years of evolution
    1664                   qsurf_hist(ig,igcm_h2o_ice) = &
    1665                      !qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0
    1666                      qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
     1542                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
    16671543
    16681544                  ! if ice has gone -ve, set to zero
    16691545                  if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
    16701546                     qsurf_hist(ig,igcm_h2o_ice) = 0.0
    1671                      !qsurf_hist(ig,igcm_h2o_vap) = 0.0
    16721547                  endif
    16731548
    … …  
    16751550                  if(ice_min(ig).lt.0.01)then
    16761551                     qsurf_hist(ig,igcm_h2o_ice) = 0.0
    1677                      !qsurf_hist(ig,igcm_h2o_vap) = 0.0
    16781552                  endif
    16791553
    … …  
    16811555
    16821556               ! enforce ice conservation
    1683                ice_tot=0.0
    1684                do ig = 1, ngrid
    1685                   ice_tot = ice_tot + qsurf_hist(ig,igcm_h2o_ice)*area(ig)
    1686                enddo
    1687                do ig = 1, ngrid
    1688                   qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice)*(icesrf/ice_tot)
    1689                enddo
     1557               ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
     1558               qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
    16901559
    16911560            endif
    … …  
    18031672             call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
    18041673             call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
     1674             call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat)
    18051675             call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
    18061676          endif
    … …  
    18181688
    18191689!     Output aerosols
    1820         if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0)   &
    1821            call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3, &
    1822            reffrad(1,1,iaero_co2))
    1823         if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0)   &
    1824            call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3, &
    1825            reffrad(1,1,iaero_h2o))
    1826         if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0)   &
    1827            call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol', &
    1828            'um kg m^-2',2,reffcol(1,iaero_co2))
    1829         if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0)   &
    1830            call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol', &
    1831            'um kg m^-2',2,reffcol(1,iaero_h2o))
     1690        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
     1691          call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
     1692        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
     1693          call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
     1694        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
     1695          call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
     1696        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
     1697          call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
    18321698
    18331699!     Output tracers
    … …  
    18431709
    18441710          if(watercond.or.CLFvarying)then
    1845 !             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
    1846 !             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
    1847 !             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
     1711             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
     1712             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
     1713             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
    18481714             call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
    18491715          endif
    … …  
    18771743!      output spectrum
    18781744      if(specOLR.and.corrk)then     
    1879                  call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
    1880                  call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
     1745         call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
     1746         call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
    18811747      endif
    18821748

  • trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90

    r726 r728  
    5151!               can in princple be anything: currently it's H2O.
    5252!
    53 !     NAERKIND  The number of radiatively active species
    54 !               (set in scatterers.h ; built when compiling with makegcm -s #)
     53!     NAERKIND  The number of radiatively active aerosol types
     54!
    5555!     NSIZEMAX  The maximum number of aerosol particle sizes
    5656!
    … …  
    7676      ! equivalent temperatures are 1/NTfac of these values
    7777      integer, parameter :: NTstar = 500
    78       integer, parameter :: NTstop = 9000 ! new default for all non hot Jupiter runs
     78      integer, parameter :: NTstop = 15000 ! new default for all non hot Jupiter runs
    7979      real*8, parameter :: NTfac = 1.0D+1 
    8080      !integer, parameter :: NTstar = 1000

  • trunk/LMDZ.GENERIC/libf/phystd/rain.F90

    r650 r728  
    22
    33
    4   use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2
    5 
     4! to use  'getin'
     5  use ioipsl_getincom
     6  use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
     7  use radii_mod, only: h2o_cloudrad
    68  implicit none
    79
    … …  
    1517!     -------
    1618!     Adapted from the LMDTERRE code by R. Wordsworth (2009)
     19!     Added rain vaporization in case of T>Tsat
    1720!     Original author Z. X. Li (1993)
    1821!     
    … …  
    5154      PARAMETER (seuil_neb=0.001)
    5255
    53       REAL,PARAMETER :: cl=2.0e-4                          ! Precipitation threshold
    54       REAL,PARAMETER :: ct=1./1800.                        ! Precipitation rate
     56      INTEGER,save :: precip_scheme      ! id number for precipitaion scheme
     57!     for simple scheme  (precip_scheme=1)
     58      REAL,SAVE :: rainthreshold                ! Precipitation threshold in simple scheme
     59!     for sundquist scheme  (precip_scheme=2-3)
     60      REAL,SAVE :: cloud_sat                    ! Precipitation threshold in non simple scheme
     61      REAL,SAVE :: precip_timescale             ! Precipitation timescale
     62!     for Boucher scheme  (precip_scheme=4)
     63      REAL,SAVE :: Cboucher ! Precipitation constant in Boucher 95 scheme
     64      REAL,PARAMETER :: Kboucher=1.19E8
     65      REAL,SAVE :: c1
    5566
    5667      INTEGER ninter
    5768      PARAMETER (ninter=5)
    58 
    59       logical simple                    ! Use very simple Emanuel scheme?
    60       parameter(simple=.true.)
    6169
    6270      logical evap_prec                 ! Does the rain evaporate?
    … …  
    6876      real lconvert
    6977
    70 !     for precipitation evaporation (old scheme)
    71       real eeff1
    72       real eeff2
    73 !      parameter (eeff1=0.95)
    74 !      parameter (eeff2=0.98)
    75       parameter (eeff1=0.5)
    76       parameter (eeff2=0.8)
    77 
    7878!     Local variables
    7979      INTEGER i, k, n
    80       REAL zqs(ngridmx,nlayermx), zdelta, zcor
     80      REAL zqs(ngridmx,nlayermx),Tsat(ngridmx,nlayermx), zdelta, zcor
    8181      REAL zrfl(ngridmx), zrfln(ngridmx), zqev, zqevt
    8282
    … …  
    8585      REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx)
    8686
    87       real ttemp, ptemp
     87      real reffh2oliq(ngridmx,nlayermx),reffh2oice(ngridmx,nlayermx)
     88 
     89      real ttemp, ptemp, psat_tmp
    8890      real tnext(ngridmx,nlayermx)
    8991
    … …  
    115117         PRINT*, 'in rain.F, ninter=', ninter
    116118         PRINT*, 'in rain.F, evap_prec=', evap_prec
     119
     120         write(*,*) "Precipitation scheme to use?"
     121         precip_scheme=1 ! default value
     122         call getin("precip_scheme",precip_scheme)
     123         write(*,*) " precip_scheme = ",precip_scheme
     124
     125         if (precip_scheme.eq.1) then
     126            write(*,*) "rainthreshold in simple scheme?"
     127            rainthreshold=0. ! default value
     128            call getin("rainthreshold",rainthreshold)
     129            write(*,*) " rainthreshold = ",rainthreshold
     130
     131         else if (precip_scheme.eq.2.or.precip_scheme.eq.3) then
     132            write(*,*) "cloud water saturation level in non simple scheme?"
     133            cloud_sat=2.6e-4   ! default value
     134            call getin("cloud_sat",cloud_sat)
     135            write(*,*) " cloud_sat = ",cloud_sat
     136            write(*,*) "precipitation timescale in non simple scheme?"
     137            precip_timescale=3600.  ! default value
     138            call getin("precip_timescale",precip_timescale)
     139            write(*,*) " precip_timescale = ",precip_timescale
     140
     141         else if (precip_scheme.eq.4) then
     142            write(*,*) "multiplicative constant in Boucher 95 precip scheme"
     143            Cboucher=1.   ! default value
     144            call getin("Cboucher",Cboucher)
     145            write(*,*) " Cboucher = ",Cboucher 
     146            c1=1.00*1.097/rhowater*Cboucher*Kboucher 
     147
     148         endif
    117149
    118150         firstcall = .false.
    … …  
    158190            ttemp = zt(i,k)
    159191            ptemp = pplay(i,k)
    160             call watersat(ttemp,ptemp,zqs(i,k))
     192!            call watersat(ttemp,ptemp,zqs(i,k))
     193            call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
     194            call Tsat_water(ptemp,Tsat(i,k))
    161195         ENDDO
    162196      ENDDO
    … …  
    180214               IF (zrfl(i) .GT.0.) THEN
    181215
    182                   zqev     = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! BC modif here
    183                   zqevt    = 2.0e-5*(1.0-q(i,k)/zqs(i,k))    &
    184                        *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here
    185                   zqevt    = MAX (zqevt, 0.0)
    186                   zqev     = MIN (zqev, zqevt)
    187                   zqev     = MAX (zqev, 0.0)
    188                   zrfln(i) = zrfl(i) - zqev
    189                   zrfln(i) = max(zrfln(i),0.0)
    190 
    191                   !zqev     = MAX (0.0, (zqs(i,k)-q(i,k))*eeff1 )
    192                   !zqevt    = (zrfl(i)/l2c(i,k))*eeff2
    193                   !zqev     = MIN (zqev, zqevt)
    194                   !zrfln(i) = zrfl(i) - zqev*l2c(i,k)
    195                   !zrfln(i)  = zrfl(i) - 1.5e-5*(1.0-q(i,k)/zqs(i,k))*sqrt(zrfl(i))
    196                   !zrfln(i)  = min(zrfln(i),zrfl(i))
    197                   ! this is what is actually written in the manual
    198 
    199                   d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
    200                   !d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
    201                   d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
    202  
    203                   zrfl(i)  = zrfln(i)
     216                  if(zt(i,k).gt.Tsat(i,k))then
     217!                    print*,'in rain',i,k,zrfl(i),zt(i,k),Tsat(i,k)
     218                     zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*l2c(i,k)/RLVTT,zrfl(i))
     219                     zrfl(i)=MAX(zrfl(i)-zqev,0.)
     220                     d_q(i,k)=zqev/l2c(i,k)
     221                     d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
     222!                    print*,zqev,zrfl(i),d_q(i,k),d_t(i,k)
     223                  else
     224                     zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! there is a bug here
     225                     zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k))    & !default was 2.e-5
     226                        *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here
     227                     zqevt = MAX (zqevt, 0.0)
     228                     zqev  = MIN (zqev, zqevt)
     229                     zqev  = MAX (zqev, 0.0)
     230                     zrfln(i)= zrfl(i) - zqev
     231                     zrfln(i)= max(zrfln(i),0.0)
     232
     233                     d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
     234                     !d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
     235                     d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
     236                     zrfl(i)  = zrfln(i)
     237                  end if
     238                     
     239
    204240               ENDIF
    205241            ENDDO
    … …  
    211247
    212248
    213          if(simple)then
     249         if(precip_scheme.eq.1)then
    214250
    215251            DO i = 1, ngridmx
    … …  
    234270            ENDDO
    235271
    236          else
    237 
     272         elseif (precip_scheme.ge.2) then
     273         
    238274            DO i = 1, ngridmx
    239275               IF (rneb(i,k).GT.0.0) THEN
    … …  
    248284            ENDDO
    249285
     286           SELECT CASE(precip_scheme)
     287 !precip scheme from Sundquist 78
     288           CASE(2)
     289
    250290            DO n = 1, ninter
    251291               DO i = 1, ngridmx
    252292                  IF (rneb(i,k).GT.0.0) THEN
    253                      zchau(i) = (ct*ptimestep/FLOAT(ninter)) * zoliq(i)      &
    254                           * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) * zfrac(i)
    255 
    256                      ! this is the ONLY place where zneb, ct and cl are used
     293                     ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
     294
     295                     zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale)) * zoliq(i)      &
     296                          * (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)**2)) * zfrac(i)
    257297                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
    258298                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
    … …  
    268308               ENDDO
    269309            ENDDO           
     310
     311 !precip scheme modified from Sundquist 78 (in q**3)
     312           CASE(3)         
     313           
     314            DO n = 1, ninter
     315               DO i = 1, ngridmx
     316                  IF (rneb(i,k).GT.0.0) THEN
     317                     ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
     318
     319                     zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale*cloud_sat**2)) * (zoliq(i)/zneb(i))**3 
     320                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
     321                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
     322                          *fallv(zrhol(i)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
     323!                          * 0.1 * (1.0-zfrac(i))
     324                     ztot(i)  = zchau(i) + zfroi(i)
     325
     326                     IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
     327                     ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
     328                     zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
     329
     330                  ENDIF
     331               ENDDO
     332            ENDDO           
     333
     334 !precip scheme modified from Boucher 95
     335           CASE(4)
     336
     337            !recalculate liquid water particle radii
     338            call h2o_cloudrad(ql,reffh2oliq,reffh2oice)
     339
     340            DO n = 1, ninter
     341               DO i = 1, ngridmx
     342                  IF (rneb(i,k).GT.0.0) THEN
     343                     ! this is the ONLY place where zneb and c1 are used
     344
     345                     zchau(i) = ptimestep/FLOAT(ninter) *c1* zrho(i) &
     346                                    *(zoliq(i)/zneb(i))**2*reffh2oliq(i,k)*zneb(i)* zfrac(i)
     347                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
     348                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
     349                          *fallv(zrhol(i)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
     350!                          * 0.1 * (1.0-zfrac(i))
     351                     ztot(i)  = zchau(i) + zfroi(i)
     352
     353                     IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
     354                     ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
     355                     zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
     356
     357                  ENDIF
     358               ENDDO
     359            ENDDO           
     360
     361           END SELECT ! precip_scheme
    270362
    271363            ! Change in cloud density and surface H2O values
    … …  
    277369            ENDDO
    278370
    279          endif ! if simple
     371
     372         endif ! if precip_scheme=1
    280373
    281374 9999 continue

  • trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F

    r726 r728  
    2626!     B. Charnay
    2727!     A. Spiga
     28!     J. Leconte (2012)
    2829!
    2930!==================================================================

  • trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90

    r253 r728  
    11      subroutine su_watercycle
    22
     3! to use  'getin'
     4      use ioipsl_getincom
    35      use watercommon_h
    46      implicit none
    57#include "comcstfi.h"
     8#include "callkeys.h"
    69
    710
    … …  
    1013!     Purpose
    1114!     -------
    12 !     Set up relevant constants and parameters for the water cycle
     15!     Set up relevant constants and parameters for the water cycle, and water cloud properties
    1316!
    1417!     Authors
    1518!     -------
    1619!     Robin Wordsworth (2010)
     20!     Jeremy Leconte (2012)
    1721!
    1822!==================================================================

  • trunk/LMDZ.GENERIC/libf/phystd/totalcloudfrac.F90

    r253 r728  
    1       subroutine totalcloudfrac(rneb,totalrneb)
     1      subroutine totalcloudfrac(rneb,totalrneb,pplev,pq,tau)
    22
     3      use watercommon_h
    34      implicit none
    45
    … …  
    2223#include "comgeomfi.h"
    2324#include "comdiurn.h"
     25#include "callkeys.h"
     26
    2427
    2528      real rneb(ngridmx,nlayermx)             ! cloud fraction     
     29      real pplev(ngridmx,nlayermx+1)
     30      real pq(ngridmx,nlayermx,nqmx)
     31      real tau(ngridmx,nlayermx)
     32
    2633      real totalrneb(ngridmx)                 ! total cloud fraction
    27 
    28       integer recovery
    29       parameter(recovery=1)
     34      real, dimension(nlayermx+1) :: masse
     35      integer, parameter          :: recovery=7
     36      integer ltau_max
     37      real massetot
    3038
    3139! hypothesis behind recovery. value:
    … …  
    3341! 2 = maximal recovery
    3442! 3 = minimal recovery
    35 
    36 
     43! 4 = fixed recovery
     44! 5 = recovery on the thicker layer
    3745!     Local variables
    3846      integer ig, l
    39       real clear
     47      real clear,tau_min
     48      real, parameter ::  tau_c=0.1 !threshold of optical depth for the calculation of total cloud fraction
     49      real rneb2(nlayermx)
     50
    4051
    4152      do ig=1,ngridmx
    … …  
    5162         elseif (recovery.eq.2) then
    5263            totalrneb(ig) = rneb(ig,1)
    53             do l=2,nlayermx       
     64            do l=2,14 !nlayermx       
    5465               totalrneb(ig) = max(rneb(ig,l),totalrneb(ig))
    5566            enddo
    … …  
    6071               totalrneb(ig) = min(rneb(ig,l),totalrneb(ig))
    6172            enddo
    62          endif                  ! (recovery=)
    6373         
     74         elseif (recovery.eq.4) then
     75            totalrneb(ig) = CLFfixval
     76
     77         elseif (recovery.eq.5) then
     78            totalrneb(ig) = rneb(ig,1)           
     79            do l=1,nlayermx
     80               masse(l)=pq(ig,l,igcm_h2o_ice)*(pplev(ig,l)-pplev(ig,l+1))
     81            enddo
     82            ltau_max=maxloc(masse,dim=1)
     83            totalrneb(ig) = rneb(ig,ltau_max)
     84
     85         elseif (recovery.eq.6) then
     86            totalrneb(ig) = 0.           
     87            do l=1,nlayermx
     88               masse(l)=pq(ig,l,igcm_h2o_ice)*(pplev(ig,l)-pplev(ig,l+1))
     89               masse(l)=max(masse(l),0.)
     90            enddo
     91            massetot=sum(masse,dim=1)
     92            do l=1,nlayermx
     93               totalrneb(ig) = totalrneb(ig)+rneb(ig,l)*masse(l)/massetot
     94            enddo
     95
     96         elseif (recovery.eq.7) then
     97
     98            rneb2(:)=rneb(ig,1:nlayermx)
     99            tau_min=MIN(tau_c,MAXVAL(tau(ig,1:nlayermx))/2.)
     100            do l=1,nlayermx
     101               if(tau(ig,l)<tau_min) rneb2(l)=0.       
     102            enddo
     103            totalrneb(ig)=maxval(rneb2(1:nlayermx))
     104
     105         endif                  ! (recovery=)   
     106
    64107         totalrneb(ig) = min(1.,totalrneb(ig))
    65108         totalrneb(ig) = max(0.,totalrneb(ig))

  • trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90

    r726 r728  
    55          pdufi,pdvfi,pdtfi,pdqfi,pfluxsrf,            &
    66          Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
    7           pdqdif,pdqsdif,lastcall)
    8 
    9       use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
     7          pdqdif,pdqevap,pdqsdif,lastcall)
     8
     9      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
    1010      use radcommon_h, only : sigma
    1111
    … …  
    112112      INTEGER iq
    113113      REAL zq(ngridmx,nlayermx,nqmx)
     114      REAL zqnoevap(ngridmx,nlayermx) !special for water case to compute where evaporated water goes.
     115      REAL pdqevap(ngridmx,nlayermx) !special for water case to compute where evaporated water goes.
     116      REAL zdmassevap(ngridmx)
    114117      REAL rho(ngridmx)         ! near-surface air density
    115       REAL qsat(ngridmx)
     118      REAL qsat(ngridmx),psat(ngridmx)
    116119      DATA firstcall/.true./
    117120      REAL kmixmin
    … …  
    119122!     Variables added for implicit latent heat inclusion
    120123!     --------------------------------------------------
    121       real dqsat(ngridmx), qsat_temp1, qsat_temp2
    122 
    123       integer ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
    124       save ivap, iliq, iliq_surf,iice_surf
     124      real dqsat(ngridmx), qsat_temp1, qsat_temp2,psat_temp
     125
     126      integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
    125127
    126128      real, parameter :: karman=0.4
    … …  
    203205            ENDDO
    204206         ENDDO
     207         if (water) then
     208            DO ilev=1,nlay
     209               DO ig=1,ngrid
     210                  zqnoevap(ig,ilev)=pq(ig,ilev,ivap) + pdqfi(ig,ilev,ivap)*ptimestep
     211               ENDDO
     212            ENDDO
     213         Endif
    205214      end if
    206215
    … …  
    246255      end if
    247256
    248 !JL12 change zkv at the surface by zcdv to calculate the surface momentum properly
     257!JL12 change zkv at the surface by zcdv to calculate the surface momentum flux properly
    249258      DO ig=1,ngrid
    250259         zkv(ig,1)=zcdv(ig)
    … …  
    253262 
    254263!JL12 calculate the flux coefficients (tables multiplied elements by elements)
    255       zfluxv=zkv(:,1:nlay)*zb0
     264      zfluxv(1:ngridmx,1:nlayermx)=zkv(1:ngridmx,1:nlayermx)*zb0(1:ngridmx,1:nlayermx)
    256265     
    257266!-----------------------------------------------------------------------
    … …  
    344353         zcdv(ig) = zcdv_true(ig)*sqrt(zu2)
    345354         zcdh(ig) = zcdh_true(ig)*sqrt(zu2)
    346          zkh(ig,1)=zcdh(ig)
     355         zkh(ig,1)= zcdh(ig)
    347356      ENDDO
    348357
    349358!     JL12 calculate the flux coefficients (tables multiplied elements by elements)
    350359!     ---------------------------------------------------------------
    351       zfluxq=zkh(:,1:nlay)*zb0 !JL12 we save zfluxq which doesn't need the Exner factor
    352       zfluxt=zfluxq*zExner
     360      zfluxq(1:ngridmx,1:nlayermx)=zkh(1:ngridmx,1:nlayermx)*zb0(1:ngridmx,1:nlayermx) !JL12 we save zfluxq which doesn't need the Exner factor
     361      zfluxt(1:ngridmx,1:nlayermx)=zfluxq(1:ngridmx,1:nlayermx)*zExner(1:ngridmx,1:nlayermx)
    353362
    354363
    … …  
    502511
    503512                ! Calculate the value of qsat at the surface (water)
    504                 call watersat(ptsrf(ig),pplev(ig,1),qsat(ig))
    505                 call watersat(ptsrf(ig)-0.0001,pplev(ig,1),qsat_temp1)
    506                 call watersat(ptsrf(ig)+0.0001,pplev(ig,1),qsat_temp2)
     513                call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
     514                call Psat_water(ptsrf(ig)-0.0001,pplev(ig,1),psat_temp,qsat_temp1)
     515                call Psat_water(ptsrf(ig)+0.0001,pplev(ig,1),psat_temp,qsat_temp2)
    507516                dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
    508517                ! calculate dQsat / dT by finite differences
    509518                ! we cannot use the updated temperature value yet...
    510  
     519!               if(qsat2(ig).gt.1.) then
     520!                  qsat(ig)=1.
     521!                  dqsat(ig)=0.
     522!               else
     523!                  qsat(ig)=qsat2(ig)
     524!               endif
    511525               enddo
     526
     527               call writediagfi(ngrid,'qsatused','saturation mixing ratio at surface','',2,qsat)
     528               call writediagfi(ngrid,'psat','saturation pressure at surface','',2,psat)
    512529
    513530! coefficients for q
    … …  
    632649           endif                ! if (water et iq=ivap)
    633650        end do                  ! of do iq=1,nq
    634       endif                     ! traceur
     651
     652        if (water) then  ! special case where we recompute water mixing without any evaporation.
     653                         !    The difference with the first calculation then tells us where evaporated water has gone
     654
     655            DO ig=1,ngrid
     656               z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
     657               zcq(ig,nlay)=zmass(ig,nlay)*zqnoevap(ig,nlay)*z1(ig)
     658               zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
     659            ENDDO
     660           
     661            DO ilay=nlay-1,2,-1
     662               DO ig=1,ngrid
     663                  z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
     664                  zcq(ig,ilay)=(zmass(ig,ilay)*zqnoevap(ig,ilay)+zfluxq(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig)
     665                  zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
     666               ENDDO
     667            ENDDO
     668
     669            DO ig=1,ngrid
     670               z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
     671               zcq(ig,1)=(zmass(ig,1)*zqnoevap(ig,1)+zfluxq(ig,2)*zcq(ig,2))*z1(ig)
     672            enddo
     673
     674!     Starting upward calculations for simple tracer mixing (e.g., dust)
     675            do ig=1,ngrid
     676               zqnoevap(ig,1)=zcq(ig,1)
     677            end do
     678
     679            do ilay=2,nlay
     680               do ig=1,ngrid
     681                  zqnoevap(ig,ilay)=zcq(ig,ilay)+zdq(ig,ilay)*zqnoevap(ig,ilay-1)
     682               end do
     683            end do
     684
     685         endif               ! if water
     686       
     687       
     688      endif                     ! tracer
    635689
    636690
    … …  
    659713            enddo
    660714         enddo
     715         if (water) then
     716            do ilev = 1, nlay
     717               do ig=1,ngrid
     718                  pdqevap(ig,ilev)=(zq(ig,ilev,ivap)-zqnoevap(ig,ilev))/ptimestep
     719               enddo
     720            enddo
     721            do ig=1,ngrid
     722               zdmassevap(ig)=SUM(pdqevap(ig,:)*zmass(ig,:))*ptimestep
     723            end do         
     724         endif
    661725      endif
    662726
    663727      if(water)then
    664       call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
     728         call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
     729         if (tracer) then
     730            call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap)
     731            call writediagfi(ngrid,'h2oevap','evaporated water vapor mass','kg/m2',2,zdmassevap)
     732         endif
    665733      endif
    666734

  • trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90

    r650 r728  
    1       module watercommon_h
     1module watercommon_h
    22
    33      implicit none
    … …  
    1111      real, parameter :: RLVTT = 2.257E+6 ! Latent heat of vaporization (J kg-1)
    1212      real, parameter :: RLSTT = 2.257E+6 ! 2.591E+6 in reality ! Latent heat of sublimation (J kg-1)
    13       !real, parameter :: RLVTT = 0.0
    14       !real, parameter :: RLSTT = 0.0
    1513
    1614      real, parameter :: RLFTT = 3.334E+5 ! Latent heat of fusion (J kg-1) ! entails an energy sink but better description of albedo
    17       !real, parameter :: RLFTT = 0.0
    1815      real, parameter :: rhowater = 1.0E+3 ! mass of water (kg/m^3)
     16      real, parameter :: rhowaterice = 9.2E+2 ! mass of water (kg/m^3)
    1917      real, parameter :: capcal_h2o_liq = 4181.3 ! specific heat capacity of liquid water J/kg/K
    2018      real, parameter :: mx_eau_sol = 150 ! mass of water (kg/m^2)
    2119
    22 !     This was the old Martian latent heat version:
    23 !     lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
     20      real, save :: epsi, RCPD, RCPV, RV, RVTMP2
     21     
     22      contains
    2423
    25       real, save :: epsi, RCPD, RCPV, RV, RVTMP2
     24     
     25!==================================================================
     26      subroutine Psat_water(T,p,psat,qsat)
    2627
    27       end module watercommon_h
     28         implicit none
     29
     30!==================================================================
     31!     Purpose
     32!     -------
     33!     Compute the saturation vapor pressure and mass mixing ratio at saturation (kg/kg)
     34!     for a given pressure (Pa) and temperature (K)
     35!     Based on the Tetens formula from L.Li physical parametrization manual
     36!
     37!     Authors
     38!     -------
     39!     Jeremy Leconte (2012)
     40!
     41!==================================================================
     42
     43!        input
     44         real, intent(in) :: T, p
     45 
     46!        output
     47         real psat,qsat
     48
     49! JL12 variables for tetens formula
     50         real,parameter :: Pref_solid_liquid=611.14
     51         real,parameter :: Trefvaporization=35.86
     52         real,parameter :: Trefsublimation=7.66
     53         real,parameter :: r3vaporization=17.269
     54         real,parameter :: r3sublimation=21.875
     55
     56! checked vs. old watersat data 14/05/2012 by JL.
     57
     58         if (T.lt.T_h2O_ice_liq) then ! solid / vapour
     59            psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation))
     60         else                 ! liquid / vapour
     61            psat = Pref_solid_liquid*Exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization))
     62         endif
     63         if(psat.gt.p) then
     64            qsat=1.
     65         else
     66            qsat=epsi*psat/(p-(1.-epsi)*psat)
     67         endif
     68         return
     69      end subroutine Psat_water
     70
     71
     72
     73
     74!==================================================================
     75      subroutine Lcpdqsat_water(T,p,psat,qsat,dqsat)
     76
     77         implicit none
     78
     79!==================================================================
     80!     Purpose
     81!     -------
     82!     Compute L/cp*d (q_sat)/d T
     83!     for a given temperature (K)!
     84!     Based on the Tetens formula from L.Li physical parametrization manual
     85!
     86!     Authors
     87!     -------
     88!     Jeremy Leconte (2012)
     89!
     90!==================================================================
     91
     92!        input
     93         real T, p, psat, qsat
     94 
     95!        output
     96         real dqsat
     97
     98! JL12 variables for tetens formula
     99         real,parameter :: Pref_solid_liquid=611.14
     100         real,parameter :: Trefvaporization=35.86
     101         real,parameter :: Trefsublimation=7.66
     102         real,parameter :: r3vaporization=17.269
     103         real,parameter :: r3sublimation=21.875
     104
     105         real :: dummy
     106
     107         if (psat.gt.p) then
     108            dqsat=0.
     109            return
     110         endif
     111
     112         if (T.lt.T_h2O_ice_liq) then ! solid / vapour
     113            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2
     114         else                 ! liquid / vapour
     115            dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2
     116         endif
     117
     118         dqsat=RLVTT/RCPD*qsat**2*(p/(epsi*psat))*dummy
     119         return
     120      end subroutine Lcpdqsat_water
     121
     122
     123
     124
     125!==================================================================
     126      subroutine Tsat_water(p,Tsat)
     127
     128         implicit none
     129
     130!==================================================================
     131!     Purpose
     132!     -------
     133!     Compute the saturation temperature
     134!     for a given pressure (Pa)
     135!     Based on the Tetens formula from L.Li physical parametrization manual
     136!
     137!     Authors
     138!     -------
     139!     Jeremy Leconte (2012)
     140!
     141!==================================================================
     142
     143!        input
     144         real p
     145 
     146!        output
     147         real Tsat
     148
     149! JL12 variables for tetens formula
     150         real,parameter :: Pref_solid_liquid=611.14
     151         real,parameter :: Trefvaporization=35.86
     152         real,parameter :: Trefsublimation=7.66
     153         real,parameter :: r3vaporization=17.269
     154         real,parameter :: r3sublimation=21.875
     155
     156         if (p.lt.Pref_solid_liquid) then ! solid / vapour
     157            Tsat =(T_h2O_ice_liq*r3sublimation- Trefsublimation*Log(p/Pref_solid_liquid))/(r3sublimation-Log(p/Pref_solid_liquid))
     158         else                 ! liquid / vapour
     159            Tsat =(T_h2O_ice_liq*r3vaporization- Trefvaporization*Log(p/Pref_solid_liquid))/(r3vaporization-Log(p/Pref_solid_liquid))
     160         endif
     161
     162         return
     163      end subroutine Tsat_water
     164
     165
     166end module watercommon_h

  • trunk/LMDZ.GENERIC/libf/phystd/watersat.F90

    r650 r728  
    2929
    3030  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
    31      qsat = 100.0 * epsi * 10**(2.07023 - 0.00320991             &
     31     qsat = 100.0 * 10**(2.07023 - 0.00320991             &
    3232          * T - 2484.896 / T + 3.56654 * alog10(T))
    3333  else                 ! liquid / vapour
    34      qsat = 100.0 * epsi * 10**(23.8319 - 2948.964 / T - 5.028  &
     34     qsat = 100.0 * 10**(23.8319 - 2948.964 / T - 5.028  &
    3535          * alog10(T) - 29810.16 * exp( -0.0699382 * T)  &
    3636          + 25.21935 * exp(-2999.924/T))
    3737  endif
    38   qsat=qsat/p
    39 
     38!  qsat=epsi*qsat/p
     39  if(qsat.gt.p) then
     40     qsat=1.
     41  else
     42     qsat=epsi*qsat/(p-(1.-epsi)*qsat)
     43  endif
    4044  return
    4145end subroutine watersat
Note: See TracChangeset for help on using the changeset viewer.