source: trunk/LMDZ.MARS/libf/phymars/simpleclouds.F @ 3026

Last change on this file since 3026 was 2378, checked in by lrossi, 5 years ago

MARS GCM

HDO
Correction of an error in newstart for inihdo.
Other minor corrections for HDO cycle.
Transition from fractionation coefficients from Merlivat et al. 1967 to Lamb et al. 2017

LR

File size: 8.9 KB
Line 
1      subroutine simpleclouds(ngrid,nlay,ptimestep,
2     &             pplay,pzlay,pt,pdt,
3     &             pq,pdq,pdqcloud,pdtcloud,
4     &             nq,tau,rice)
5      USE updaterad
6      USE watersat_mod, ONLY: watersat
7      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
8     &                      igcm_hdo_vap, igcm_hdo_ice,
9     &                      qperemin
10      USE comcstfi_h
11      use dimradmars_mod, only: naerkind
12
13      implicit none
14c------------------------------------------------------------------
15c  This routine is used to form clouds when a parcel of the GCM is
16c    saturated. It is a simplified scheme, and there is almost no
17c    microphysics involved. When the air is saturated, water-ice
18c    clouds form on a fraction of the dust particles, specified by
19c    the constant called "ccn_factor". There is no supersaturation,
20c    and no nucleation rates computed. A more accurate scheme can
21c    be found in the routine called "improvedclouds.F".
22
23c  Modif de zq si saturation dans l'atmosphere
24c  si zq(ig,l)> zqsat(ig,l) ->    zq(ig,l)=zqsat(ig,l)
25c  Le test est effectue de bas en haut. L'eau condensee
26c    (si saturation) est remise dans la couche en dessous.
27c  L'eau condensee dans la couche du bas est deposee a la surface
28
29c  Authors: Franck Montmessin (water ice scheme)
30c           Francois Forget (changed nuclei density & outputs)
31c           Ehouarn Millour (sept.2008, tracers are now handled
32c                                   by name and not fixed index)
33c           J.-B. Madeleine (developed a single routine called
34c                            simpleclouds.F, and corrected calculations
35c                            of the typical CCN profile, Oct. 2011)
36c------------------------------------------------------------------
37#include "callkeys.h"
38
39c------------------------------------------------------------------
40c     Arguments:
41c     ---------
42c     Inputs:
43      INTEGER ngrid,nlay
44      integer nq                 ! nombre de traceurs
45      REAL ptimestep             ! pas de temps physique (s)
46      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
47      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
48      REAL pt(ngrid,nlay)        ! temperature at the middle of the
49                                 !   layers (K)
50      REAL pdt(ngrid,nlay)       ! tendance temperature des autres
51                                 !   param.
52      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
53      real pdq(ngrid,nlay,nq)    ! tendance avant condensation
54                                 !   (kg/kg.s-1)
55      REAL tau(ngrid,naerkind)   ! Column dust optical depth at each point
56
57c     Output:
58      REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
59                                 ! (r_c in montmessin_2004)
60      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
61                                   !   H2O(kg/kg.s-1)
62      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
63                                   !   a la chaleur latente
64
65c------------------------------------------------------------------
66c     Local variables:
67
68      REAL rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
69
70      INTEGER ig,l
71
72      REAL zq(ngrid,nlay,nq)    ! local value of tracers
73      REAL zq0(ngrid,nlay,nq)   ! local initial value of tracers
74      REAL zt(ngrid,nlay)       ! local value of temperature
75      REAL zqsat(ngrid,nlay)    ! saturation
76      REAL*8 dzq                      ! masse de glace echangee (kg/kg)
77      REAL lw                         !Latent heat of sublimation (J.kg-1)
78      REAL,PARAMETER :: To=273.15     ! reference temperature, T=273.15 K
79      real rdusttyp(ngrid,nlay) ! Typical dust geom. mean radius (m)
80      REAL ccntyp(ngrid,nlay)
81                                      ! Typical dust number density (#/kg)
82      REAL alpha_c(ngrid,nlay) !!MARGAUX: alpha_c as a spatial variable
83
84c     CCN reduction factor
85c      REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant
86      REAL DoH_vap(ngrid,nlay)
87      REAL DoH_ice(ngrid,nlay)
88c-----------------------------------------------------------------------
89c    1. initialisation
90c    -----------------
91
92c    On "update" la valeur de q(nq) (water vapor) et temperature.
93c    On effectue qqes calculs preliminaires sur les couches :
94
95      do l=1,nlay
96        do ig=1,ngrid
97          zq(ig,l,igcm_h2o_vap)=
98     &      pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
99          zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
100          zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
101          zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
102
103          zq(ig,l,igcm_h2o_ice)=
104     &      pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
105          zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004
106          zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
107
108          if (hdo) then
109          zq(ig,l,igcm_hdo_vap)=
110     &      pq(ig,l,igcm_hdo_vap)+pdq(ig,l,igcm_hdo_vap)*ptimestep
111          zq(ig,l,igcm_hdo_vap)=max(zq(ig,l,igcm_hdo_vap),1e-30) ! FF 12/2004
112          zq0(ig,l,igcm_hdo_vap)=zq(ig,l,igcm_hdo_vap)
113
114          zq(ig,l,igcm_hdo_ice)=
115     &      pq(ig,l,igcm_hdo_ice)+pdq(ig,l,igcm_hdo_ice)*ptimestep
116          zq(ig,l,igcm_hdo_ice)=max(zq(ig,l,igcm_hdo_ice),1e-30) ! FF 12/2004
117          zq0(ig,l,igcm_hdo_ice)=zq(ig,l,igcm_hdo_ice)
118
119          endif !hdo
120        enddo
121      enddo
122
123      pdqcloud(1:ngrid,1:nlay,1:nq)=0
124      pdtcloud(1:ngrid,1:nlay)=0
125
126      alpha_c(1:ngrid,1:nlay)=0.
127
128c     ----------------------------------------------
129c
130c     Rapport de melange a saturation dans la couche l : -------
131c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132
133      call watersat(ngrid*nlay,zt,pplay,zqsat)
134
135c     taux de condensation (kg/kg/s-1) dans les differentes couches
136c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
137
138      do l=1,nlay
139        do ig=1,ngrid
140
141          if (zq(ig,l,igcm_h2o_vap).ge.zqsat(ig,l))then  !  Condensation
142            dzq=zq(ig,l,igcm_h2o_vap)-zqsat(ig,l)
143
144          elseif(zq(ig,l,igcm_h2o_vap).lt.zqsat(ig,l))then  ! Sublimation
145            dzq=-min(zqsat(ig,l)-zq(ig,l,igcm_h2o_vap),
146     &               zq(ig,l,igcm_h2o_ice))
147          endif
148           
149c         Water Mass change
150c         ~~~~~~~~~~~~~~~~~
151          zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)+dzq
152          zq(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)-dzq
153       
154        enddo ! of do ig=1,ngrid
155      enddo ! of do l=1,nlay
156
157
158c     Tendance finale
159c     ~~~~~~~~~~~~~~~
160      do l=1, nlay
161        do ig=1,ngrid
162          pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
163     &                            -zq0(ig,l,igcm_h2o_vap))/ptimestep
164          pdqcloud(ig,l,igcm_h2o_ice) =
165     &      (zq(ig,l,igcm_h2o_ice) - zq0(ig,l,igcm_h2o_ice))/ptimestep
166
167          lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
168          pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
169
170          if (hdo) then
171            if (pdqcloud(ig,l,igcm_h2o_ice).gt.0.0) then !condens
172
173                if (hdofrac) then ! do we use fractionation?
174c               alpha_c(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
175                alpha_c(ig,l) = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
176                else
177                alpha_c(ig,l) = 1.d0
178                endif
179               
180                if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
181                  pdqcloud(ig,l,igcm_hdo_ice)=
182     &              pdqcloud(ig,l,igcm_h2o_ice)*alpha_c(ig,l)*
183     &         ( zq0(ig,l,igcm_hdo_vap)
184     &                 /zq0(ig,l,igcm_h2o_vap) )
185                else
186                   pdqcloud(ig,l,igcm_hdo_ice)= 0.0
187                endif
188
189                pdqcloud(ig,l,igcm_hdo_ice) =
190     &            min(pdqcloud(ig,l,igcm_hdo_ice),
191     &              zq0(ig,l,igcm_hdo_vap)/ptimestep)
192
193                pdqcloud(ig,l,igcm_hdo_vap)=
194     &               -pdqcloud(ig,l,igcm_hdo_ice)       
195
196          else  ! sublimation
197
198             if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
199                pdqcloud(ig,l,igcm_hdo_ice)=
200     &               pdqcloud(ig,l,igcm_h2o_ice)*
201     &      ( zq0(ig,l,igcm_hdo_ice)
202     &              /zq0(ig,l,igcm_h2o_ice) )
203             else
204                pdqcloud(ig,l,igcm_hdo_ice)= 0.
205             endif
206
207              pdqcloud(ig,l,igcm_hdo_ice) =
208     &          max(pdqcloud(ig,l,igcm_hdo_ice),
209     &            -zq0(ig,l,igcm_hdo_ice)/ptimestep)
210
211              pdqcloud(ig,l,igcm_hdo_vap)=
212     &             -pdqcloud(ig,l,igcm_hdo_ice)       
213
214            endif ! condensation/sublimation
215
216          endif ! hdo
217
218        enddo ! of do ig=1,ngrid
219      enddo ! of do l=1,nlay
220
221c     ice crystal radius
222      do l=1, nlay
223        do ig=1,ngrid
224          call updaterice_typ(zq(ig,l,igcm_h2o_ice),
225     &       tau(ig,1),pzlay(ig,l),rice(ig,l))
226        end do
227      end do
228
229c     if (hdo) then
230c           CALL WRITEDIAGFI(ngrid,'alpha_c',
231c    &                       'alpha_c',
232c    &                       ' ',3,alpha_c)
233c     endif !hdo
234c------------------------------------------------------------------
235      return
236      end
Note: See TracBrowser for help on using the repository browser.