Ignore:
Timestamp:
May 4, 2012, 6:35:58 PM (13 years ago)
Author:
jleconte
Message:
  • Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
  • Minor name changes in watercommon
  • Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

  • Added consistency check in inifis
  • Moved 1d water initialization from physiqu to rcm1d
Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
14 edited

Legend:

Unmodified
Added
Removed

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

    r622 r650  
    172172
    173173!     are particle radii fixed?
    174       logical radfixed
     174      logical, save :: radfixed
    175175      real maxrad, minrad
    176176
     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     
    177184      real CBRT
    178185      external CBRT
    … …  
    183190      REAL*8 muvarrad(L_LEVELS)
    184191
    185       radfixed=.false.
    186 
    187192!=======================================================================
    188193!     Initialization on first call
    … …  
    195200      qsiaer(:,:,:)=0.0
    196201      giaer(:,:,:) =0.0
     202      radfixed=.false.
    197203
    198204
    199205      if(firstcall) then
     206         if(kastprof)then
     207            radfixed=.true.
     208         endif 
    200209
    201210         call system('rm -f surf_vals_long.out')
    … …  
    305314      enddo
    306315
    307       if(kastprof)then
    308          radfixed=.true.
    309       endif
     316
    310317
    311318      if(radfixed)then
    … …  
    320327               do ig=1,ngrid
    321328                  !reffrad(ig,l,2) = 2.e-5 ! H2O ice
    322                   reffrad(ig,l,2) = 5.e-6 ! H2O ice
     329                   reffrad(ig,l,2) = rad_chaud ! H2O ice
     330
     331                   zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
     332                   zfice = MIN(MAX(zfice,0.0),1.0)
     333                   reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice
     334                   nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice
    323335               enddo
    324336            enddo

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

    r596 r650  
    1111     &   , Nmix_co2, Nmix_h2o                                           &
    1212     &   , dusttau                                                      &
    13      &   , nonideal                                                     &
    1413     &   , meanOLR                                                      &
    1514     &   , specOLR                                                      &

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

    r253 r650  
    6161         ! ignoring qevap term creates huge difference when qevap large!!!
    6262
    63          zdelta = MAX(0.,SIGN(1.,To-tevap(ig,l))) ! what is this?
     63         zdelta = MAX(0.,SIGN(1.,T_h2O_ice_liq-tevap(ig,l))) ! what is this?
    6464                                                  ! for division between water / liquid
    6565         dqevap(ig,l) = MAX(0.0,qevap(ig,l,igcm_h2o_ice))/dtime

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

    r368 r650  
    33     albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
    44
    5   use watercommon_h, only: To, RLFTT, rhowater, mx_eau_sol
     5  use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
    66
    77  implicit none
    … …  
    189189            hice(ig)    = 0.
    190190
    191             if(twater(ig) .lt. To)then
    192                E=min((To+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick)
     191            if(twater(ig) .lt. T_h2O_ice_liq)then
     192               E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick)
    193193               hice(ig)        = E/(RLFTT*rhowater)
    194194               hice(ig)        = max(hice(ig),0.0)
    … …  
    221221
    222222!     melt the snow
    223             if(ztsurf(ig).gt.To)then
     223            if(ztsurf(ig).gt.T_h2O_ice_liq)then
    224224               if(zqsurf(ig,iice).gt.1.0e-8)then
    225225
    226                   a     = (ztsurf(ig)-To)*pcapcal(ig)/RLFTT
     226                  a     = (ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
    227227                  b     = zqsurf(ig,iice)
    228228                  fsnoi = min(a,b)
    … …  
    240240               if(zqsurf(ig,iliq).gt.1.0e-8)then
    241241
    242                   a     = -(ztsurf(ig)-To)*pcapcal(ig)/RLFTT
     242                  a     = -(ztsurf(ig)-T_h2O_ice_liq)*pcapcal(ig)/RLFTT
    243243                  b     = zqsurf(ig,iliq)
    244244                 

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

    r622 r650  
    219219         call getin("calladj",calladj)
    220220         write(*,*) " calladj = ",calladj
    221        
    222 
    223 ! Test of incompatibility
    224          if (enertest.and.nonideal) then
    225             print*,'Energy conservation calculations currently
    226      &           assume ideal gas!'
    227             call abort
    228          endif
    229221
    230222         write(*,*) "call CO2 condensation ?"
    … …  
    232224         call getin("co2cond",co2cond)
    233225         write(*,*) " co2cond = ",co2cond
    234 
     226! Test of incompatibility
     227         if (co2cond.and.(.not.tracer)) then
     228            print*,'We need a CO2 ice tracer to condense CO2'
     229            call abort
     230         endif   
     231         
    235232         write(*,*) "CO2 supersaturation level ?"
    236233         co2supsat=1.0 ! default value

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

    r253 r650  
    11subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb)
    22
    3   use watercommon_h, only: To, RLVTT, RCPD
     3  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD
    44
    55  implicit none
    … …  
    124124         DO i = 1, ngridmx
    125125            v_cptt(i,k) = RCPD * local_t(i,k)
    126          ENDDO
    127       ENDDO
    128 
    129       DO k = 1, nlayermx
    130          DO i = 1, ngridmx
    131             v_cptt(i,k) = RCPD * local_t(i,k)
    132126            v_t = local_t(i,k)
    133127            v_p = pplay(i,k)
    … …  
    142136!         DO i = 1, ngridmx
    143137!            v_t = local_t(i,k)
    144 !            IF (v_t.LT.To) THEN
     138!            IF (v_t.LT.T_h2O_ice_liq) THEN
    145139!               print*,'RHs=',q(i,k) / v_qs(i,k)
    146140!            ELSE

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

    r622 r650  
    88 
    99      use radinc_h, only : naerkind,L_NSPECTI,L_NSPECTV
    10       use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
     10      use watercommon_h, only : RLVTT
    1111      use gases_h
    1212      use radcommon_h, only: sigma
    … …  
    370370      real vdifcncons(ngridmx)
    371371      real cadjncons(ngridmx)
    372       real qzero1D
    373       save qzero1D
    374372
    375373!      double precision qsurf_hist(ngridmx,nqmx)
    … …  
    544542            ! initialise variables for water cycle
    545543
    546             qzero1D=0.0
    547 
    548544            if(enertest)then
    549545               watertest = .true.
    550             endif
    551 
    552             if(ngrid.eq.1)then
    553                qzero1D               = 0.0
    554                qsurf(1,igcm_h2o_vap) = qzero1D
    555                do l=1, nlayer
    556                   pq(1,l,igcm_h2o_vap)=0.0
    557                   pq(1,l,igcm_h2o_ice)=0.0
    558                enddo
    559546            endif
    560547

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

    r622 r650  
    22
    33
    4   use watercommon_h, only: To, RLVTT, RCPD, RCPV, RV, RVTMP2
     4  use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2
    55
    66  implicit none
    … …  
    8282
    8383      REAL zoliq(ngridmx)
    84       REAL ztglace
    8584      REAL zdz(ngridmx),zrho(ngridmx),ztot(ngridmx), zrhol(ngridmx)
    8685      REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx)
    … …  
    140139      ENDDO
    141140      ENDDO
    142 
    143 !     Determine the cold clouds by their temperature
    144       ztglace = To - 15.0
    145141
    146142!     Initialise the outputs
    … …  
    219215            DO i = 1, ngridmx
    220216               ttemp = zt(i,k)
    221                IF (ttemp .ge. To) THEN
     217               IF (ttemp .ge. T_h2O_ice_liq) THEN
    222218                  lconvert=rainthreshold
    223219               ELSEIF (ttemp .gt. t_crit) THEN
    … …  
    245241                  zrho(i)  = pplay(i,k) / ( zt(i,k) * R )
    246242                  zdz(i)   = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
    247                   zfrac(i) = (zt(i,k)-ztglace) / (To-ztglace)
     243                  zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
    248244                  zfrac(i) = MAX(zfrac(i), 0.0)
    249245                  zfrac(i) = MIN(zfrac(i), 1.0)
    … …  
    291287            call abort
    292288         endif
    293          IF (t(i,1) .LT. To) THEN
     289         IF (t(i,1) .LT. T_h2O_ice_liq) THEN
    294290            dqssnow(i) = zrfl(i)
    295291            dqsrain(i) = 0.0

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

    r601 r650  
    7373!      REAL co2ice           ! co2ice layer (kg.m-2) !not used anymore
    7474      integer :: i_co2_ice=0  ! tracer index of co2 ice
    75       integer :: i_h2o_ice=0  ! tracer index of co2 ice
     75      integer :: i_h2o_ice=0  ! tracer index of h2o ice
     76      integer :: i_h2o_vap=0  ! tracer index of h2o vapor
    7677      REAL emis             ! surface layer
    7778      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
    … …  
    111112      REAL cloudfrac(ngridmx,nlayermx)
    112113      REAL hice(ngridmx),totcloudfrac(ngridmx)
     114      REAL qzero1D   !initial water amount on the ground
    113115
    114116c=======================================================================
    … …  
    228230         i_co2_ice=0
    229231         i_h2o_ice=0
     232         i_h2o_vap=0
    230233         do iq=1,nq
    231234           if (tnom(iq)=="co2_ice") then
    … …  
    233236           elseif (tnom(iq)=="h2o_ice") then
    234237             i_h2o_ice=iq
     238           elseif (tnom(iq)=="h2o_vap") then
     239             i_h2o_vap=iq
    235240           endif
    236241         enddo
    … …  
    275280      area(1)=1.E+0
    276281      aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
    277      !!! GEOPOTENTIAL. useless in 1D because control bu surface pressure
     282     !!! GEOPOTENTIAL. useless in 1D because control by surface pressure
    278283      phisfi(1)=0.E+0
    279284     !!! LATITUDE. can be set.
    … …  
    493498        qsurf(iq) = 0.
    494499      ENDDO
     500
     501      if (water) then
     502         qzero1D               = 0.0
     503         qsurf(i_h2o_vap) = qzero1D
     504      endif
    495505
    496506c   Initialisation pour prendre en compte les vents en 1-D

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

    r600 r650  
    77          pdqdif,pdqsdif,lastcall)
    88
    9       use watercommon_h, only : RLVTT, To, RCPD, mx_eau_sol
     9      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
    1010      use radcommon_h, only : sigma
    1111
    … …  
    609609                     else !dqsdif_total(ig).ge.0
    610610                        !If water vapour is condensing, we must decide whether it forms ice or liquid.
    611                         if(ztsrf(ig).gt.To)then
     611                        if(ztsrf(ig).gt.T_h2O_ice_liq)then
    612612                           pdqsdif(ig,iice_surf)=0.0
    613613                           pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep

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

    r600 r650  
    77     &     pdqdif,pdqsdif,lastcall)
    88
    9       use watercommon_h, only : RLVTT, To, RCPD, mx_eau_sol
     9      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
    1010      use radcommon_h, only : sigma
    1111
    … …  
    644644!     If water vapour is condensing, we must decide whether it forms ice or liquid.
    645645                 if(dqsdif_total(ig).gt.0)then ! a bug was here!
    646                     if(ztsrf2(ig).gt.To)then
     646                    if(ztsrf2(ig).gt.T_h2O_ice_liq)then
    647647                       pdqsdif(ig,iice)=0.0
    648648                       pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
    … …  
    704704               Wdiff=Wtot/ptimestep+pdqsdif(ig,ivap)+pdqsdif(ig,iice)
    705705
    706                if(ztsrf2(ig).gt.To)then
     706               if(ztsrf2(ig).gt.T_h2O_ice_liq)then
    707707                  pdqsdif(ig,ivap)=pdqsdif(ig,ivap)-Wdiff
    708708               else

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

    r253 r650  
    44
    55      real, parameter :: T_coup = 234.0
    6       real, parameter :: To = 273.16
     6      real, parameter :: T_h2O_ice_liq = 273.16
     7      real, parameter :: T_h2O_ice_clouds = T_h2O_ice_liq-15.
    78      real, parameter :: mH2O = 18.01528   
    89
    910      ! benjamin additions
    10       real, parameter :: RLVTT = 2.257E+6 ! Latent heat of sublimation (J kg-1)
     11      real, parameter :: RLVTT = 2.257E+6 ! Latent heat of vaporization (J kg-1)
    1112      real, parameter :: RLSTT = 2.257E+6 ! 2.591E+6 in reality ! Latent heat of sublimation (J kg-1)
    1213      !real, parameter :: RLVTT = 0.0
    1314      !real, parameter :: RLSTT = 0.0
    1415
    15       real, parameter :: RLFTT = 3.334E+5 ! Latent heat of fusion (J kg-1)
     16      real, parameter :: RLFTT = 3.334E+5 ! Latent heat of fusion (J kg-1) ! entails an energy sink but better description of albedo
    1617      !real, parameter :: RLFTT = 0.0
    1718      real, parameter :: rhowater = 1.0E+3 ! mass of water (kg/m^3)
     19      real, parameter :: capcal_h2o_liq = 4181.3 ! specific heat capacity of liquid water J/kg/K
    1820      real, parameter :: mx_eau_sol = 150 ! mass of water (kg/m^2)
    1921
    … …  
    2123!     lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
    2224
    23       real epsi, RCPD, RCPV, RV, RVTMP2
    24       save epsi, RCPD, RCPV, RV, RVTMP2
     25      real, save :: epsi, RCPD, RCPV, RV, RVTMP2
    2526
    2627      end module watercommon_h

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

    r253 r650  
    11subroutine watersat(T,p,qsat)
    22
    3   use watercommon_h, only: To, epsi
     3  use watercommon_h, only: T_h2O_ice_liq, epsi
    44  implicit none
    55
    … …  
    2828! x by epsi converts to mass mixing ratio
    2929
    30   if (T.lt.To) then ! solid / vapour
     30  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
    3131     qsat = 100.0 * epsi * 10**(2.07023 - 0.00320991             &
    3232          * T - 2484.896 / T + 3.56654 * alog10(T))

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

    r135 r650  
    11subroutine watersat_grad(T,qsat,dqsat)
    22
    3   use watercommon_h, only: T_coup, RLVTT, RCPD
     3  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD,T_coup
    44  implicit none
    55
    … …  
    77!     Purpose
    88!     -------
    9 !     Compute the water mass mixing ratio at saturation (kg/kg)
    10 !     for a given pressure (Pa) and temperature (K)
     9!     Compute the L/cp*d (q_sat)/d T
     10!     for a given temperature (K)
    1111!
    1212!     Authors
    … …  
    2222  real dqsat
    2323
    24   if (T.lt.T_coup) then ! solid / vapour
     24!  if (T.lt.T_coup) then ! solid / vapour !why use T_coup?????????? JL12
     25  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
    2526     dqsat = RLVTT/RCPD*qsat*(3.56654/T             &
    2627          +2484.896*LOG(10.)/T**2                   &
Note: See TracChangeset for help on using the changeset viewer.