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

Last change on this file since 2613 was 2073, checked in by aslmd, 6 years ago

in largescale enforce type conversion as is done in the other instance of Lcpdqsat_water. corrects an error with picky compilers (e.g. gfortran)
File size: 6.5 KB
Line 
1      subroutine largescale(ngrid,nlayer,nq,ptimestep, pplev, pplay,    &
2                    pt, pq, pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
3
4
5      use ioipsl_getin_p_mod, only: getin_p
6      use watercommon_h, only : RLVTT, RCPD, RVTMP2,  &
7          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water,Lcpdqsat_water
8      USE tracer_h
9      IMPLICIT none
10
11!==================================================================
12!     
13!     Purpose
14!     -------
15!     Calculates large-scale (stratiform) H2O condensation.
16!     
17!     Authors
18!     -------
19!     Adapted from the LMDTERRE code by R. Wordsworth (2009)
20!     Original author Z. X. Li (1993)
21!     
22!==================================================================
23
24      INTEGER ngrid,nlayer,nq
25
26!     Arguments
27      REAL ptimestep                 ! intervalle du temps (s)
28      REAL pplev(ngrid,nlayer+1) ! pression a inter-couche
29      REAL pplay(ngrid,nlayer)   ! pression au milieu de couche
30      REAL pt(ngrid,nlayer)      ! temperature (K)
31      REAL pq(ngrid,nlayer,nq) ! tracer mixing ratio (kg/kg)
32      REAL pdt(ngrid,nlayer)     ! physical temperature tenedency (K/s)
33      REAL pdq(ngrid,nlayer,nq)! physical tracer tenedency (K/s)
34      REAL pdtlsc(ngrid,nlayer)  ! incrementation de la temperature (K)
35      REAL pdqvaplsc(ngrid,nlayer) ! incrementation de la vapeur d'eau
36      REAL pdqliqlsc(ngrid,nlayer) ! incrementation de l'eau liquide
37      REAL rneb(ngrid,nlayer)    ! fraction nuageuse
38
39
40!     Options du programme
41      REAL, SAVE :: ratqs   ! determine largeur de la distribution de vapeur
42!$OMP THREADPRIVATE(ratqs)
43
44!     Variables locales
45      REAL CBRT
46      EXTERNAL CBRT
47      INTEGER i, k , nn
48      INTEGER,PARAMETER :: nitermax=5000
49      DOUBLE PRECISION,PARAMETER :: alpha=.1,qthreshold=1.d-8
50      ! JL13: if "careful, T<Tmin in psat water" appears often, you may want to stabilise the model by
51      !                   decreasing alpha and increasing nitermax accordingly
52      DOUBLE PRECISION zq(ngrid)
53      DOUBLE PRECISION zcond(ngrid),zcond_iter
54      DOUBLE PRECISION zdelq(ngrid)
55      DOUBLE PRECISION zqs(ngrid)
56      real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs
57     
58! evaporation calculations
59      REAL dqevap(ngrid,nlayer),dtevap(ngrid,nlayer)     
60      REAL qevap(ngrid,nlayer,nq)
61      REAL tevap(ngrid,nlayer)
62
63      DOUBLE PRECISION zx_q(ngrid)
64      LOGICAL,SAVE :: firstcall=.true.
65!$OMP THREADPRIVATE(firstcall)
66
67
68      IF (firstcall) THEN
69
70         write(*,*) "value for ratqs? "
71         ratqs=0.2 ! default value
72         call getin_p("ratqs",ratqs)
73         write(*,*) " ratqs = ",ratqs
74
75         firstcall = .false.
76      ENDIF
77
78!     GCM -----> subroutine variables, initialisation of outputs
79
80      pdtlsc(1:ngrid,1:nlayer)  = 0.0
81      pdqvaplsc(1:ngrid,1:nlayer)  = 0.0
82      pdqliqlsc(1:ngrid,1:nlayer) = 0.0
83      rneb(1:ngrid,1:nlayer) = 0.0
84      Lcp=RLVTT/RCPD
85
86
87      ! Evaporate cloud water/ice
88      call evap(ngrid,nlayer,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
89      ! note: we use qevap but not tevap in largescale/moistadj
90            ! otherwise is a big mess
91
92
93!  Boucle verticale (du haut vers le bas)
94   DO k = nlayer, 1, -1
95
96      zt(1:ngrid)=pt(1:ngrid,k)+(pdt(1:ngrid,k)+dtevap(1:ngrid,k))*ptimestep
97      zq(1:ngrid)=qevap(1:ngrid,k,igcm_h2o_vap) !liquid water is included in qevap
98
99!     Calculer la vapeur d'eau saturante et
100!     determiner la condensation partielle
101      DO i = 1, ngrid
102
103         local_p=pplay(i,k)
104         if(zt(i).le.15.) then
105            print*,'in lsc',i,k,zt(i)
106!           zt(i)=15.   ! check too low temperatures
107         endif
108         call Psat_water(zt(i),local_p,psat_tmp,zqs_temp)
109         zqs(i)=zqs_temp
110 
111         zdelq(i) = MAX(MIN(ratqs * zq(i),1.-zq(i)),1.d-12)
112         rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
113#ifdef MESOSCALE
114         if (rneb(i,k).lt.0.01) then  !no clouds MESO
115#else
116         if (rneb(i,k).lt.0.) then  !no clouds
117#endif
118
119            rneb(i,k)=0.
120            zcond(i)=0.
121
122         else if ((rneb(i,k).gt.0.99).or.(ratqs.lt.1.e-6)) then    !complete cloud cover, we start without evaporating
123            rneb(i,k)=1.
124            zt(i)=pt(i,k)+pdt(i,k)*ptimestep
125            zx_q(i) = pq(i,k,igcm_h2o_vap)+pdq(i,k,igcm_h2o_vap)*ptimestep
126            dqevap(i,k)=0.
127!           iterative process to stabilize the scheme when large water amounts JL12
128            zcond(i) = 0.0d0
129            Do nn=1,nitermax 
130               call Psat_water(zt(i),local_p,psat_tmp,zqs_temp)
131               zqs(i)=zqs_temp
132               call Lcpdqsat_water(zt(i),local_p,psat_tmp,zqs_temp,zdqs,dlnpsat_tmp)
133               zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs)     
134                  !zcond can be negative here
135               zx_q(i) = zx_q(i) - zcond_iter
136               zcond(i) = zcond(i) + zcond_iter
137               zt(i) = zt(i) + zcond_iter*Lcp
138               if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit
139!              if (ABS(zcond_iter/alpha).lt.qthreshold) exit
140               if (nn.eq.nitermax) print*,'itermax in largescale'
141            End do ! niter
142            zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_h2o_ice)+pdq(i,k,igcm_h2o_ice)*ptimestep))
143
144         else   !standard case     
145#ifdef MESOSCALE
146            rneb(i,k)=1. !LES/MESO case
147            zx_q(i) = (zq(i)+zqs(i))/2.0d0 ! LES
148#else
149            zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0d0 !water vapor in cloudy sky
150#endif
151!           iterative process to stabilize the scheme when large water amounts JL12
152            zcond(i) = 0.0d0
153            Do nn=1,nitermax
154               ! use zqs_temp and not zqs(i) to force type conversion
155               ! -- might not be a good solution, actually
156               ! but this is compliant with "complete cloud cover" case above
157               call Lcpdqsat_water(zt(i),local_p,psat_tmp,zqs_temp,zdqs,dlnpsat_tmp)
158               zcond_iter = MAX(0.0d0,alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs))         
159                  !zcond always postive! cannot evaporate clouds!
160                  !this is why we must reevaporate before largescale
161               zx_q(i) = zx_q(i) - zcond_iter
162               zcond(i) = zcond(i) + zcond_iter
163               if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit
164!              if (ABS(zcond_iter/alpha).lt.qthreshold) exit
165               zt(i) = zt(i) + zcond_iter*Lcp*rneb(i,k)
166               call Psat_water(zt(i),local_p,psat_tmp,zqs_temp)
167               zqs(i)=zqs_temp
168               if (nn.eq.nitermax) print*,'itermax in largescale'
169            End do ! niter
170
171         Endif
172
173         zcond(i) = zcond(i)*rneb(i,k)/ptimestep ! JL12
174
175      ENDDO
176
177!     Tendances de t et q
178         pdqvaplsc(1:ngrid,k)  = dqevap(1:ngrid,k) - zcond(1:ngrid)
179         pdqliqlsc(1:ngrid,k) = - pdqvaplsc(1:ngrid,k)
180         pdtlsc(1:ngrid,k)  = pdqliqlsc(1:ngrid,k)*Lcp
181
182   Enddo ! k= nlayer, 1, -1
183 
184
185      end
Note: See TracBrowser for help on using the repository browser.