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

Last change on this file since 2316 was 2312, checked in by lrossi, 6 years ago

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