source: trunk/LMDZ.GENERIC/libf/phystd/largescale.F90 @ 735

Last change on this file since 735 was 728, checked in by jleconte, 13 years ago

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.

File size: 4.8 KB
Line 
1      subroutine largescale(ptimestep, pplev, pplay, pt, pq,   &
2                        pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
3
4      use watercommon_h, only : RLVTT, RCPD, RVTMP2,  &
5          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water,Lcpdqsat_water
6      IMPLICIT none
7
8!==================================================================
9!     
10!     Purpose
11!     -------
12!     Calculates large-scale (stratiform) H2O condensation.
13!     
14!     Authors
15!     -------
16!     Adapted from the LMDTERRE code by R. Wordsworth (2009)
17!     Original author Z. X. Li (1993)
18!     
19!==================================================================
20
21#include "dimensions.h"
22#include "dimphys.h"
23#include "comcstfi.h"
24
25#include "fisice.h"
26#include "callkeys.h"
27#include "tracer.h"
28
29
30
31!     Arguments
32      REAL ptimestep                 ! intervalle du temps (s)
33      REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche
34      REAL pplay(ngridmx,nlayermx)   ! pression au milieu de couche
35      REAL pt(ngridmx,nlayermx)      ! temperature (K)
36      real pq(ngridmx,nlayermx,nqmx) ! tracer mixing ratio (kg/kg)
37      REAL pdt(ngridmx,nlayermx)     ! physical temperature tenedency (K/s)
38      REAL pdq(ngridmx,nlayermx,nqmx)! physical tracer tenedency (K/s)
39      REAL pdtlsc(ngridmx,nlayermx)  ! incrementation de la temperature (K)
40      REAL pdqvaplsc(ngridmx,nlayermx) ! incrementation de la vapeur d'eau
41      REAL pdqliqlsc(ngridmx,nlayermx) ! incrementation de l'eau liquide
42      REAL rneb(ngridmx,nlayermx)    ! fraction nuageuse
43
44
45!     Options du programme
46      REAL ratqs   ! determine largeur de la distribution de vapeur
47      PARAMETER (ratqs=0.2)
48
49!     Variables locales
50      REAL CBRT
51      EXTERNAL CBRT
52      INTEGER i, k , nn
53      INTEGER,PARAMETER :: niter=4
54      REAL zt(ngridmx), zq(ngridmx)
55      REAL zcond(ngridmx),zcond_iter
56      REAL zdelq(ngridmx)
57      REAL zqs(ngridmx), zdqs(ngridmx)
58      REAL psat_tmp
59     
60! evaporation calculations
61      REAL dqevap(ngridmx,nlayermx),dtevap(ngridmx,nlayermx)     
62      REAL qevap(ngridmx,nlayermx,nqmx)
63      REAL tevap(ngridmx,nlayermx)
64
65      REAL zcor(ngridmx), zdelta(ngridmx), zcvm5(ngridmx)
66      REAL zx_q(ngridmx)
67      REAL Nmix_local,zfice
68
69!     GCM -----> subroutine variables, initialisation of outputs
70
71      pdtlsc(1:ngridmx,1:nlayermx)  = 0.0
72      pdqvaplsc(1:ngridmx,1:nlayermx)  = 0.0
73      pdqliqlsc(1:ngridmx,1:nlayermx) = 0.0
74      rneb(1:ngridmx,1:nlayermx) = 0.0
75
76
77      ! Evaporate cloud water/ice
78      call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
79      ! note: we use qevap but not tevap in largescale/moistadj
80            ! otherwise is a big mess
81
82
83!  Boucle verticale (du haut vers le bas)
84   DO k = nlayermx, 1, -1
85
86      zt(1:ngridmx)=pt(1:ngridmx,k)+(pdt(1:ngridmx,k)+dtevap(1:ngridmx,k))*ptimestep
87      zq(1:ngridmx)=qevap(1:ngridmx,k,igcm_h2o_vap) !liquid water is included in qevap
88
89!     Calculer la vapeur d'eau saturante et
90!     determiner la condensation partielle
91      DO i = 1, ngridmx
92
93         if(zt(i).le.15.) zt(i)=15.   ! check too low temperatures
94!         call watersat(zt(i),pplay(i,k),zqs(i))
95         call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
96 
97         zdelq(i) = ratqs * zq(i)
98         if(zq(i)+zdelq(i).lt.0.999) then
99            rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
100            zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0 !water vapor in cloudy sky
101            if (rneb(i,k) .LE. 0.0) zx_q(i) = 0.0
102            if (rneb(i,k) .GE. 1.0) zx_q(i) = zq(i)
103            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
104         else
105            if(zq(i).gt.zqs(i)) then
106               rneb(i,k)=1.
107               zx_q(i)=zq(i)
108            else
109               rneb(i,k)=0.
110               zx_q(i)=zqs(i)  !no condensation needed
111            Endif
112         Endif
113
114!        iterative process to stabilize the scheme when large water amounts JL12
115         zcond(i) = 0.0
116         Do nn=1,niter 
117!            call watersat_grad(zt(i),zqs(i),zdqs(i))
118            call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
119            zcond_iter = MAX(0.0,(zx_q(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i)))         
120               !zcond always postive! cannot evaporate clouds!
121               !this is why we must reevaporate before largescale
122            zx_q(i) = zx_q(i) - zcond_iter
123            zcond(i) = zcond(i) + zcond_iter
124            zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
125            call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
126         End do ! niter
127
128
129         zcond(i) = zcond(i)/ptimestep ! added by RDW
130      ENDDO
131
132!     Tendances de t et q
133         pdqvaplsc(1:ngridmx,k)  = dqevap(1:ngridmx,k) - zcond(1:ngridmx)
134         pdqliqlsc(1:ngridmx,k) = - pdqvaplsc(1:ngridmx,k)
135         pdtlsc(1:ngridmx,k)  = pdqliqlsc(1:ngridmx,k)*RLVTT/RCPD
136
137   Enddo ! k= nlayermx, 1, -1
138   
139      !print*,'qsat=',zqs
140      !print*,'q=',q
141      !print*,'dq=',pdqvaplsc*ptimestep
142      !print*,'dT in LS=',pdtlsc*ptimestep
143
144      !print*,'rice=',rice
145      !print*,'rneb=',rneb
146
147      return
148      end
Note: See TracBrowser for help on using the repository browser.