source: trunk/LMDZ.MARS/libf/phymars/watercloud.F @ 1448

Last change on this file since 1448 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

  • Optimized computations in paramfoto_compact (twice less dlog10 calculations)
  • Checked consistency before/after modification in debug mode
  • Checked performance is not impacted (same as before)
File size: 15.3 KB
Line 
1       SUBROUTINE watercloud(ngrid,nlay,ptimestep,
2     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
3     &                pq,pdq,pdqcloud,pdtcloud,
4     &                nq,tau,tauscaling,rdust,rice,nuice,
5     &                rsedcloud,rhocloud)
6! to use  'getin'
7      USE ioipsl_getincom
8      USE updaterad
9      USE comcstfi_h
10      use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
11     &                      igcm_dust_mass, igcm_dust_number,
12     &                      igcm_ccn_mass, igcm_ccn_number,
13     &                      rho_dust, nuice_sed, nuice_ref
14      use dimradmars_mod, only: naerkind
15      IMPLICIT NONE
16
17
18c=======================================================================
19c  Water-ice cloud formation
20c 
21c  Includes two different schemes:
22c    - A simplified scheme (see simpleclouds.F)
23c    - An improved microphysical scheme (see improvedclouds.F)
24c
25c  There is a time loop specific to cloud formation
26c  due to timescales smaller than the GCM integration timestep.
27c
28c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
29c           J.-B. Madeleine, Thomas Navarro
30c
31c  2004 - 2012
32c=======================================================================
33
34c-----------------------------------------------------------------------
35c   declarations:
36c   -------------
37
38#include "callkeys.h"
39
40c   Inputs:
41c   ------
42
43      INTEGER ngrid,nlay
44      INTEGER nq                 ! nombre de traceurs
45      REAL ptimestep             ! pas de temps physique (s)
46      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
47      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
48      REAL pdpsrf(ngrid)         ! tendence surf pressure
49      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
50      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
51      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
52
53      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
54      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
55
56      REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
57      REAL tauscaling(ngrid)   ! Convertion factor for dust amount
58      real rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
59
60c   Outputs:
61c   -------
62
63      real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
64      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
65                                   ! a la chaleur latente
66
67      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
68                               ! (r_c in montmessin_2004)
69      REAL nuice(ngrid,nlay)   ! Estimated effective variance
70                               !   of the size distribution
71      real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
72      real rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
73
74c   local:
75c   ------
76     
77      ! for ice radius computation
78      REAL Mo,No
79      REAl ccntyp
80     
81      ! for time loop
82      INTEGER microstep  ! time subsampling step variable
83      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
84      SAVE imicro
85      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
86      SAVE microtimestep
87     
88      ! tendency given by clouds (inside the micro loop)
89      REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
90      REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
91
92      ! global tendency (clouds+physics)
93      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
94      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
95
96
97      INTEGER iq,ig,l
98      LOGICAL,SAVE :: firstcall=.true.
99
100c    ** un petit test de coherence
101c       --------------------------
102
103      IF (firstcall) THEN
104         
105        if (nq.gt.nqmx) then
106           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
107           write(*,*) 'nq=',nq,' nqmx=',nqmx
108           stop
109        endif
110         
111        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
112        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
113               
114        write(*,*) "time subsampling for microphysic ?"
115#ifdef MESOSCALE
116        imicro = 2
117#else
118        imicro = 30
119#endif
120        call getin("imicro",imicro)
121        write(*,*)"imicro = ",imicro
122       
123        microtimestep = ptimestep/real(imicro)
124        write(*,*)"Physical timestep is",ptimestep
125        write(*,*)"Microphysics timestep is",microtimestep
126
127        firstcall=.false.
128      ENDIF ! of IF (firstcall)
129     
130c-----Initialization
131      subpdq(1:ngrid,1:nlay,1:nq) = 0
132      subpdt(1:ngrid,1:nlay)      = 0
133     
134      ! default value if no ice
135      rhocloud(1:ngrid,1:nlay) = rho_dust
136
137
138
139c------------------------------------------------------------------
140c Time subsampling for microphysics
141c------------------------------------------------------------------
142      DO microstep=1,imicro
143     
144c-------------------------------------------------------------------
145c   1.  Tendencies:
146c------------------
147
148
149c------ Temperature tendency subpdt
150        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
151        ! If imicro=1 subpdt is the same as pdt
152        DO l=1,nlay
153          DO ig=1,ngrid
154             subpdt(ig,l) = subpdt(ig,l)
155     &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
156          ENDDO
157        ENDDO
158c------ Tracers tendencies subpdq
159c------ At each micro timestep we add pdq in order to have a stepped entry
160        IF (microphys) THEN
161          DO l=1,nlay
162            DO ig=1,ngrid
163              subpdq(ig,l,igcm_dust_mass) =
164     &            subpdq(ig,l,igcm_dust_mass)
165     &          + pdq(ig,l,igcm_dust_mass)
166              subpdq(ig,l,igcm_dust_number) =
167     &            subpdq(ig,l,igcm_dust_number)
168     &          + pdq(ig,l,igcm_dust_number)
169              subpdq(ig,l,igcm_ccn_mass) =
170     &            subpdq(ig,l,igcm_ccn_mass)
171     &          + pdq(ig,l,igcm_ccn_mass)
172              subpdq(ig,l,igcm_ccn_number) =
173     &            subpdq(ig,l,igcm_ccn_number)
174     &          + pdq(ig,l,igcm_ccn_number)
175            ENDDO
176          ENDDO
177        ENDIF
178        DO l=1,nlay
179          DO ig=1,ngrid
180            subpdq(ig,l,igcm_h2o_ice) =
181     &          subpdq(ig,l,igcm_h2o_ice)
182     &        + pdq(ig,l,igcm_h2o_ice)
183            subpdq(ig,l,igcm_h2o_vap) =
184     &          subpdq(ig,l,igcm_h2o_vap)
185     &        + pdq(ig,l,igcm_h2o_vap)
186          ENDDO
187        ENDDO
188       
189       
190c-------------------------------------------------------------------
191c   2.  Main call to the different cloud schemes:
192c------------------------------------------------
193        IF (microphys) THEN
194           CALL improvedclouds(ngrid,nlay,microtimestep,
195     &             pplay,pt,subpdt,
196     &             pq,subpdq,subpdqcloud,subpdtcloud,
197     &             nq,tauscaling)
198
199        ELSE
200           CALL simpleclouds(ngrid,nlay,microtimestep,
201     &             pplay,pzlay,pt,subpdt,
202     &             pq,subpdq,subpdqcloud,subpdtcloud,
203     &             nq,tau,rice)
204        ENDIF
205
206
207c-------------------------------------------------------------------
208c   3.  Updating tendencies after cloud scheme:
209c-----------------------------------------------
210
211        IF (microphys) THEN
212          DO l=1,nlay
213            DO ig=1,ngrid
214              subpdq(ig,l,igcm_dust_mass) =
215     &            subpdq(ig,l,igcm_dust_mass)
216     &          + subpdqcloud(ig,l,igcm_dust_mass)
217              subpdq(ig,l,igcm_dust_number) =
218     &            subpdq(ig,l,igcm_dust_number)
219     &          + subpdqcloud(ig,l,igcm_dust_number)
220              subpdq(ig,l,igcm_ccn_mass) =
221     &            subpdq(ig,l,igcm_ccn_mass)
222     &          + subpdqcloud(ig,l,igcm_ccn_mass)
223              subpdq(ig,l,igcm_ccn_number) =
224     &            subpdq(ig,l,igcm_ccn_number)
225     &          + subpdqcloud(ig,l,igcm_ccn_number)
226            ENDDO
227          ENDDO
228        ENDIF
229        DO l=1,nlay
230          DO ig=1,ngrid
231            subpdq(ig,l,igcm_h2o_ice) =
232     &          subpdq(ig,l,igcm_h2o_ice)
233     &        + subpdqcloud(ig,l,igcm_h2o_ice)
234            subpdq(ig,l,igcm_h2o_vap) =
235     &          subpdq(ig,l,igcm_h2o_vap)
236     &        + subpdqcloud(ig,l,igcm_h2o_vap)
237          ENDDO
238        ENDDO
239       
240        IF (activice) THEN
241          DO l=1,nlay
242            DO ig=1,ngrid
243              subpdt(ig,l) =
244     &            subpdt(ig,l) + subpdtcloud(ig,l)
245            ENDDO
246          ENDDO
247        ENDIF
248     
249 
250      ENDDO ! of DO microstep=1,imicro
251     
252c-------------------------------------------------------------------
253c   6.  Compute final tendencies after time loop:
254c------------------------------------------------
255
256c------ Temperature tendency pdtcloud
257       DO l=1,nlay
258         DO ig=1,ngrid
259             pdtcloud(ig,l) =
260     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
261          ENDDO
262       ENDDO
263       
264c------ Tracers tendencies pdqcloud
265       DO l=1,nlay
266         DO ig=1,ngrid
267            pdqcloud(ig,l,igcm_h2o_ice) =
268     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro)
269     &       - pdq(ig,l,igcm_h2o_ice)
270            pdqcloud(ig,l,igcm_h2o_vap) =
271     &        subpdq(ig,l,igcm_h2o_vap)/real(imicro)
272     &       - pdq(ig,l,igcm_h2o_vap)
273         ENDDO
274       ENDDO
275       
276       IF(microphys) THEN
277        DO l=1,nlay
278         DO ig=1,ngrid
279            pdqcloud(ig,l,igcm_ccn_mass) =
280     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro)
281     &       - pdq(ig,l,igcm_ccn_mass)
282            pdqcloud(ig,l,igcm_ccn_number) =
283     &        subpdq(ig,l,igcm_ccn_number)/real(imicro)
284     &       - pdq(ig,l,igcm_ccn_number)
285         ENDDO
286        ENDDO
287       ENDIF
288       
289       IF(scavenging) THEN
290        DO l=1,nlay
291         DO ig=1,ngrid
292            pdqcloud(ig,l,igcm_dust_mass) =
293     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
294     &       - pdq(ig,l,igcm_dust_mass)
295            pdqcloud(ig,l,igcm_dust_number) =
296     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
297     &       - pdq(ig,l,igcm_dust_number)
298         ENDDO
299        ENDDO
300       ENDIF
301
302c------- Due to stepped entry, other processes tendencies can add up to negative values
303c------- Therefore, enforce positive values and conserve mass
304
305
306       IF(microphys) THEN
307        DO l=1,nlay
308         DO ig=1,ngrid
309          IF ((pq(ig,l,igcm_ccn_number) +
310     &      ptimestep* (pdq(ig,l,igcm_ccn_number) +
311     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
312     &   .or. (pq(ig,l,igcm_ccn_mass) +
313     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) +
314     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
315         pdqcloud(ig,l,igcm_ccn_number) =
316     &     - pq(ig,l,igcm_ccn_number)/ptimestep
317     &     - pdq(ig,l,igcm_ccn_number) + 1.
318         pdqcloud(ig,l,igcm_dust_number) = 
319     &     -pdqcloud(ig,l,igcm_ccn_number)
320         pdqcloud(ig,l,igcm_ccn_mass) =
321     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
322     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
323         pdqcloud(ig,l,igcm_dust_mass) =
324     &     -pdqcloud(ig,l,igcm_ccn_mass)
325          ENDIF
326         ENDDO
327        ENDDO
328       ENDIF
329
330       IF(scavenging) THEN
331        DO l=1,nlay
332         DO ig=1,ngrid
333          IF ((pq(ig,l,igcm_dust_number) +
334     &      ptimestep* (pdq(ig,l,igcm_dust_number) +
335     &        pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
336     &   .or. (pq(ig,l,igcm_dust_mass) +
337     &      ptimestep* (pdq(ig,l,igcm_dust_mass) +
338     &        pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
339         pdqcloud(ig,l,igcm_dust_number) =
340     &     - pq(ig,l,igcm_dust_number)/ptimestep
341     &     - pdq(ig,l,igcm_dust_number) + 1.
342         pdqcloud(ig,l,igcm_ccn_number) = 
343     &     -pdqcloud(ig,l,igcm_dust_number)
344         pdqcloud(ig,l,igcm_dust_mass) =
345     &     - pq(ig,l,igcm_dust_mass)/ptimestep
346     &     - pdq(ig,l,igcm_dust_mass) + 1.e-20
347         pdqcloud(ig,l,igcm_ccn_mass) =
348     &     -pdqcloud(ig,l,igcm_dust_mass)
349          ENDIF
350         ENDDO
351        ENDDO
352       ENDIF
353
354        DO l=1,nlay
355         DO ig=1,ngrid
356          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
357     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
358     &       .le. 1.e-8) THEN
359           pdqcloud(ig,l,igcm_h2o_ice) =
360     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
361           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
362          ENDIF
363          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
364     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
365     &       .le. 1.e-8) THEN
366           pdqcloud(ig,l,igcm_h2o_vap) =
367     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
368           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
369          ENDIF
370         ENDDO
371        ENDDO
372
373
374c------Update the ice and dust particle size "rice" for output or photochemistry
375c------Only rsedcloud is used for the water cycle
376
377      IF(scavenging) THEN
378        DO l=1, nlay
379         DO ig=1,ngrid
380
381        call updaterdust(
382     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
383     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
384     &    pdqcloud(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
385     &    pq(ig,l,igcm_dust_number) +                 ! dust number
386     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
387     &    pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
388     &    rdust(ig,l))
389
390         ENDDO
391        ENDDO
392      ENDIF
393       
394       
395      IF(microphys) THEN
396       
397       DO l=1, nlay
398         DO ig=1,ngrid
399
400        call updaterice_micro(
401     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
402     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
403     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
404     &    pq(ig,l,igcm_ccn_mass) +                   ! ccn mass
405     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
406     &    pdqcloud(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
407     &    pq(ig,l,igcm_ccn_number) +                 ! ccn number
408     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
409     &    pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
410     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
411         
412         ENDDO
413       ENDDO
414       
415      ELSE ! no microphys
416       
417        DO l=1,nlay
418          DO ig=1,ngrid
419         
420        call updaterice_typ(
421     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
422     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
423     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
424     &    tau(ig,1),pzlay(ig,l),rice(ig,l))
425
426          ENDDO
427         ENDDO
428       
429       ENDIF ! of IF(microphys)
430     
431     
432     
433c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
434c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
435c     Then that should not affect the ice particle radius
436      do ig=1,ngrid
437        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
438          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
439     &    rice(ig,2)=rice(ig,3)
440          rice(ig,1)=rice(ig,2)
441        end if
442      end do
443       
444       
445       DO l=1,nlay
446         DO ig=1,ngrid
447           rsedcloud(ig,l)=max(rice(ig,l)*
448     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
449     &                    rdust(ig,l))
450!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
451         ENDDO
452       ENDDO
453       
454! used for rad. transfer calculations
455! nuice is constant because a lognormal distribution is prescribed
456      nuice(1:ngrid,1:nlay)=nuice_ref
457
458
459
460c=======================================================================
461
462      END
463 
Note: See TracBrowser for help on using the repository browser.