source: trunk/LMDZ.MARS/libf/phymars/co2cloud.F @ 1635

Last change on this file since 1635 was 1629, checked in by jaudouard, 8 years ago

Bug fix concerning riceco2 declaration and some cleaning concerning the CO2 clouds scheme

  • Property svn:executable set to *
File size: 27.4 KB
Line 
1       SUBROUTINE co2cloud(ngrid,nlay,ptimestep,
2     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
3     &                pq,pdq,pdqcloudco2,pdtcloudco2,
4     &                nq,tau,tauscaling,rdust,rice,riceco2,nuice,
5     &                rsedcloudco2,rhocloudco2,zlev,pdqs_sedco2)
6! to use  'getin'
7      use dimradmars_mod, only: naerkind
8      USE comcstfi_h
9      USE ioipsl_getincom
10      USE updaterad
11      use conc_mod, only: mmean
12      use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
13     &     igcm_dust_mass, igcm_dust_number,
14     &     igcm_dust_mass, igcm_dust_number,
15     &     igcm_ccnco2_mass, igcm_ccnco2_number,
16     &     rho_dust, nuiceco2_sed, nuiceco2_ref,
17     &     rho_ice_co2,r3n_q
18      IMPLICIT NONE
19
20
21c=======================================================================
22c CO2 clouds formation
23c
24c  There is a time loop specific to cloud formation
25c  due to timescales smaller than the GCM integration timestep.
26c  microphysics subroutine is improvedCO2clouds.F
27
28c  The co2 clouds tracers (co2_ice, ccn mass and concentration) are
29c  sedimented at each microtimestep. pdqs_sedco2 keeps track of the
30c  CO2 flux at the surface
31c
32c  Authors: 09/2016 Joachim Audouard & Constantino Listowski
33c  Adaptation of the water ice clouds scheme (with specific microphysics)
34c  of Montmessin, Navarro & al.
35c
36
37c=======================================================================
38
39c-----------------------------------------------------------------------
40c   declarations:
41c   -------------
42
43!#include "dimensions.h"
44!#include "dimphys.h"
45#include "callkeys.h"
46!#include "tracer.h"
47!#include "comgeomfi.h"
48!#include "dimradmars.h"
49! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
50!#include"scatterers.h"
51#include "microphys.h"
52
53
54c   Inputs:
55c   ------
56
57      INTEGER ngrid,nlay
58      INTEGER nq                 ! nombre de traceurs
59      REAL ptimestep            ! pas de temps physique (s)
60      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
61      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
62      REAL pdpsrf(ngrid)         ! tendence surf pressure
63      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
64      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
65      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
66      real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
67
68      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
69      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
70
71      real rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
72                                ! used for nucleation of CO2 on ice-coated ccns
73
74      REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
75      REAL tauscaling(ngrid)   ! Convertion factor for dust amount
76      real rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
77
78c   Outputs:
79c   -------
80
81      real pdqcloudco2(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
82      REAL pdtcloudco2(ngrid,nlay)    ! tendence temperature due
83                                   ! a la chaleur latente
84
85      DOUBLE PRECISION riceco2(ngrid,nlay)    ! Ice mass mean radius (m)
86                               ! (r_c in montmessin_2004)
87      REAL nuice(ngrid,nlay)   ! Estimated effective variance
88                               !   of the size distribution
89      real rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
90      real rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
91      real rhocloudco2t(ngrid,nlay)  ! Cloud density (kg.m-3)
92      real  pdqs_sedco2(ngrid) ! CO2 flux at the surface
93c   local:
94c   ------
95     
96      ! for ice radius computation
97      REAL Mo,No
98      REAl ccntyp
99     
100      ! for time loop
101      INTEGER microstep  ! time subsampling step variable
102      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
103      SAVE imicro
104      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
105      SAVE microtimestep
106     
107      ! tendency given by clouds (inside the micro loop)
108      REAL subpdqcloudco2(ngrid,nlay,nq) ! cf. pdqcloud
109      REAL subpdtcloudco2(ngrid,nlay)    ! cf. pdtcloud
110
111      ! global tendency (clouds+physics)
112      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
113      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
114      real wq(ngrid,nlay+1)  !  ! displaced tracer mass (kg.m-2) during microtimestep because sedim (?/m-2)
115
116      REAL satuco2(ngrid,nlay)  ! co2 satu ratio for output
117      REAL zqsatco2(ngrid,nlay) ! saturation co2
118
119      INTEGER iq,ig,l
120      LOGICAL,SAVE :: firstcall=.true.
121      DOUBLE PRECISION Nccnco2, Niceco2,mdustJA,ndustJA
122      DOUBLE PRECISION Qccnco2
123      real :: beta
124
125      real epaisseur (ngrid,nlay) ! Layer thickness (m)
126      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
127   
128
129      real tempo_traceur_t(ngrid,nlay)
130      real tempo_traceurs(ngrid,nlay,nq)
131      real sav_trac(ngrid,nlay,nq)
132      real pdqsed(ngrid,nlay,nq)
133c     ** un petit test de coherence
134c       --------------------------
135
136      IF (firstcall) THEN
137         
138        if (nq.gt.nqmx) then
139           write(*,*) 'stop in co2cloud (nq.gt.nqmx)!'
140           write(*,*) 'nq=',nq,' nqmx=',nqmx
141           stop
142        endif
143         
144        write(*,*) "co2cloud: igcm_co2=",igcm_co2
145        write(*,*) "            igcm_co2_ice=",igcm_co2_ice
146               
147        write(*,*) "time subsampling for microphysic ?"
148#ifdef MESOSCALE
149        imicro = 2
150#else
151        imicro = 30
152#endif
153        call getin("imicro",imicro)
154             imicro=ptimestep/50.
155
156        write(*,*)"imicro = ",imicro
157       
158        microtimestep = ptimestep/real(imicro)
159        write(*,*)"Physical timestep is",ptimestep
160        write(*,*)"CO2 Microphysics timestep is",microtimestep
161
162        firstcall=.false.
163      ENDIF                     ! of IF (firstcall)
164   
165c-----Initialization
166      beta=0.85
167      subpdq(1:ngrid,1:nlay,1:nq) = 0
168      subpdt(1:ngrid,1:nlay)      = 0
169      subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
170      subpdtcloudco2(1:ngrid,1:nlay)      = 0
171     
172
173      wq(:,:)=0
174      ! default value if no ice
175      rhocloudco2(1:ngrid,1:nlay) = rho_dust
176      rhocloudco2t(1:ngrid,1:nlay) = rho_dust
177      epaisseur(1:ngrid,1:nlay)=0
178      masse(1:ngrid,1:nlay)=0
179
180      tempo_traceur_t(1:ngrid,1:nlay)=0
181      tempo_traceurs(1:ngrid,1:nlay,1:nq)=0
182      sav_trac(1:ngrid,1:nlay,1:nq)=0
183      pdqsed(1:ngrid,1:nlay,1:nq)=0
184     
185      do  l=1,nlay
186        do ig=1, ngrid
187          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
188          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
189       
190       enddo
191      enddo
192           
193     
194   
195
196           
197
198
199     
200c-------------------------------------------------------------------
201c   1.  Tendencies:
202c------------------
203 
204
205
206c------------------------------------------------------------------
207c Time subsampling for microphysics
208c------------------------------------------------------------------
209      DO microstep=1,imicro
210c------ Temperature tendency subpdt
211        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
212        ! If imicro=1 subpdt is the same as pdt
213        DO l=1,nlay
214          DO ig=1,ngrid
215c          tempo_traceur_t(ig,l)=tempo_traceur_t(ig,l)
216c     &            + subpdtcloudco2(ig,l)
217             !write(*,*) 'T micro= ', tempo_traceur_t(ig,l)
218c             tempo_traceurs(ig,l,:)=tempo_traceurs(ig,l,:)
219c     &            +subpdqcloudco2(ig,l,:)
220
221             subpdt(ig,l) = subpdt(ig,l)
222     &            + pdt(ig,l)   ! At each micro timestep we add pdt in order to have a stepped entry
223       
224             subpdq(ig,l,igcm_dust_mass) =
225     &            subpdq(ig,l,igcm_dust_mass)
226     &            + pdq(ig,l,igcm_dust_mass)
227             
228             subpdq(ig,l,igcm_dust_number) =
229     &            subpdq(ig,l,igcm_dust_number)
230     &            + pdq(ig,l,igcm_dust_number)
231             
232             subpdq(ig,l,igcm_ccnco2_mass) =
233     &            subpdq(ig,l,igcm_ccnco2_mass)
234     &            + pdq(ig,l,igcm_ccnco2_mass)
235             
236             subpdq(ig,l,igcm_ccnco2_number) =
237     &            subpdq(ig,l,igcm_ccnco2_number)
238     &            + pdq(ig,l,igcm_ccnco2_number)
239             
240             subpdq(ig,l,igcm_co2_ice) =
241     &            subpdq(ig,l,igcm_co2_ice)
242     &            + pdq(ig,l,igcm_co2_ice)
243             subpdq(ig,l,igcm_co2) =
244     &            subpdq(ig,l,igcm_co2)
245     &            + pdq(ig,l,igcm_co2)
246
247             tempo_traceur_t(ig,l)= pt(ig,l)+subpdt(ig,l)*microtimestep
248             tempo_traceurs(ig,l,:)= pq(ig,l,:)+subpdq(ig,l,:)
249     &            *microtimestep
250             !Stepped entry for sedimentation                   
251          ENDDO
252       ENDDO
253   
254!RSEDCLOUD AND RICECO2 HERE
255       
256       DO l=1, nlay
257          DO ig=1,ngrid
258             Niceco2=tempo_traceurs(ig,l,igcm_co2_ice)
259             Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number),
260     &            1.e-30)
261             Qccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_mass),
262     &            1.e-30)
263             mdustJA= tempo_traceurs(ig,l,igcm_dust_mass)
264             ndustJA=tempo_traceurs(ig,l,igcm_dust_number)
265             if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt.
266     &            1.e-30 *tauscaling(ig))) then
267                rdust(ig,l)=1.e-10
268             else
269                rdust(ig,l)=(3./4./pi/2500.*mdustJA/ndustJA)**(1./3.)
270                rdust(ig,l)=min(rdust(ig,l),5.e-6)
271                rdust(ig,l)=max(rdust(ig,l),1.e-9)   
272             end if
273             rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2
274     &            + Qccnco2*rho_dust)
275     &            / (Niceco2 + Qccnco2)
276             riceco2(ig,l)= Niceco2*3.0/
277     &            (4.0*rho_ice_co2*pi*Nccnco2)
278     &            +rdust(ig,l)*rdust(ig,l)*rdust(ig,l)
279             riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
280             write(*,*) "in co2clouds, rice = ",riceco2(ig,l)
281             write(*,*) "in co2clouds, rho = ",rhocloudco2t(ig,l)
282
283              call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
284     &            tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
285              write(*,*) "in co2clouds, rice update = ",riceco2(ig,l)
286              write(*,*) "in co2clouds, rho update = "
287     &             ,rhocloudco2t(ig,l)
288
289              rsedcloudco2(ig,l)=max(riceco2(ig,l)*
290     &            (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed),
291     &            rdust(ig,l))
292             rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),5.e-4)
293             write(*,*) 'Rsedcloud = ',rsedcloudco2(ig,l)
294             !write(*,*) 'Rhocloudco2 = ',rhocloudco2t(ig,l)
295
296          ENDDO
297       ENDDO
298       
299!     Gravitational sedimentation
300
301!     sedimentation computed from radius computed from q in module radii_mod
302       sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
303       sav_trac(:,:,igcm_ccnco2_mass)=
304     &      tempo_traceurs(:,:,igcm_ccnco2_mass)
305       sav_trac(:,:,igcm_ccnco2_number)=
306     &      tempo_traceurs(:,:,igcm_ccnco2_number)
307       
308      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
309     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
310     &     rsedcloudco2,rhocloudco2t,
311     &     tempo_traceurs(:,:,igcm_co2_ice),wq,beta) !  3 traceurs
312     
313! sedim at the surface of co2 ice
314      do ig=1,ngrid
315         pdqs_sedco2(ig)=pdqs_sedco2(ig)+  wq(ig,1)
316      end do
317
318      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
319     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
320     &     rsedcloudco2,rhocloudco2t,
321     &     tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta)
322     
323      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
324     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
325     &     rsedcloudco2,rhocloudco2t,
326     &     tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta)
327       
328     
329      DO l = 1, nlay
330         DO ig=1,ngrid
331            pdqsed(ig,l,igcm_ccnco2_mass)=
332     &           (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
333     &           sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
334            pdqsed(ig,l,igcm_ccnco2_number)=
335     &           (tempo_traceurs(ig,l,igcm_ccnco2_number)-
336     &           sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
337            pdqsed(ig,l,igcm_co2_ice)=
338     &           (tempo_traceurs(ig,l,igcm_co2_ice)-
339     &           sav_trac(ig,l,igcm_co2_ice))/microtimestep
340         ENDDO
341      ENDDO
342            !pdqsed est la tendance due a la sedimentation
343     
344      DO l = 1, nlay
345         DO ig=1,ngrid
346            pdqsed(ig,l,igcm_ccnco2_mass)=
347     &           (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
348     &           sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
349            pdqsed(ig,l,igcm_ccnco2_number)=
350     &           (tempo_traceurs(ig,l,igcm_ccnco2_number)-
351     &           sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
352            pdqsed(ig,l,igcm_co2_ice)=
353     &           (tempo_traceurs(ig,l,igcm_co2_ice)-
354     &           sav_trac(ig,l,igcm_co2_ice))/microtimestep
355         ENDDO
356      ENDDO
357            !pdqsed est la tendance due a la sedimentation
358      DO l=1,nlay
359         DO ig=1,ngrid
360            subpdq(ig,l,igcm_ccnco2_mass) =
361     &           subpdq(ig,l,igcm_ccnco2_mass)
362     &           +pdqsed(ig,l,igcm_ccnco2_mass)
363             
364            subpdq(ig,l,igcm_ccnco2_number) =
365     &           subpdq(ig,l,igcm_ccnco2_number)
366     &           +pdqsed(ig,l,igcm_ccnco2_number)
367
368            subpdq(ig,l,igcm_co2_ice) =
369     &           subpdq(ig,l,igcm_co2_ice)
370     &           +pdqsed(ig,l,igcm_co2_ice)
371         ENDDO
372      ENDDO   
373c-------------------------------------------------------------------
374c   2.  Main call to the different cloud schemes:
375c------------------------------------------------
376        IF (microphysco2) THEN
377           CALL improvedCO2clouds(ngrid,nlay,microtimestep,
378     &             pplay,pt,subpdt,
379     &             pq,subpdq,subpdqcloudco2,subpdtcloudco2,
380     &             nq,tauscaling)
381
382        ELSE
383
384           write(*,*) ' no simpleCO2clouds procedure: STOP' ! listo
385           STOP
386
387c           CALL simpleclouds(ngrid,nlay,microtimestep,   ! for water-ice clouds
388c     &             pplay,pzlay,pt,subpdt,
389c     &             pq,subpdq,subpdqcloud,subpdtcloud,
390c     &             nq,tau,riceco2)
391        ENDIF
392       
393
394c-------------------------------------------------------------------
395c   3.  Updating tendencies after cloud scheme:
396c-----------------------------------------------
397
398c        IF (microphysco2) THEN
399          DO l=1,nlay
400            DO ig=1,ngrid
401               subpdq(ig,l,igcm_dust_mass) =
402     &              subpdq(ig,l,igcm_dust_mass)
403     &              + subpdqcloudco2(ig,l,igcm_dust_mass)
404 
405               subpdq(ig,l,igcm_dust_number) =
406     &              subpdq(ig,l,igcm_dust_number)
407     &              + subpdqcloudco2(ig,l,igcm_dust_number)
408             
409               subpdq(ig,l,igcm_ccnco2_mass) =
410     &              subpdq(ig,l,igcm_ccnco2_mass)
411     &              + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
412c     &              +pdqsed(ig,l,igcm_ccnco2_mass)
413             
414               subpdq(ig,l,igcm_ccnco2_number) =
415     &              subpdq(ig,l,igcm_ccnco2_number)
416     &              + subpdqcloudco2(ig,l,igcm_ccnco2_number)
417c     &              +pdqsed(ig,l,igcm_ccnco2_number)
418
419            subpdq(ig,l,igcm_co2_ice) =
420     &           subpdq(ig,l,igcm_co2_ice)
421     &           + subpdqcloudco2(ig,l,igcm_co2_ice)
422c     &              +pdqsed(ig,l,igcm_co2_ice)
423           
424            subpdq(ig,l,igcm_co2) =
425     &           subpdq(ig,l,igcm_co2)
426     &           + subpdqcloudco2(ig,l,igcm_co2)
427         ENDDO
428        ENDDO
429       
430       
431!ici
432!      call WRITEdiagfi(ngrid,"co2cloud000","co2 traceur","kg/kg",1,
433!     &      pq(1,:,igcm_co2_ice) + ptimestep*
434!     &      ( subpdq(1,:,igcm_co2_ice)))
435     
436       
437        IF (activice) THEN
438          DO l=1,nlay
439            DO ig=1,ngrid
440              subpdt(ig,l) =
441     &            subpdt(ig,l) + subpdtcloudco2(ig,l)
442            ENDDO
443          ENDDO
444        ENDIF
445 
446 
447      ENDDO ! of DO microstep=1,imicro
448     
449c-------------------------------------------------------------------
450c   6.  Compute final tendencies after time loop:
451c------------------------------------------------
452c CO2 flux at surface (kg.m-2.s-1)
453      do ig=1,ngrid
454         pdqs_sedco2(ig)=pdqs_sedco2(ig)/ptimestep
455      enddo
456
457c------ Temperature tendency pdtcloud
458       DO l=1,nlay
459         DO ig=1,ngrid
460             pdtcloudco2(ig,l) =
461     &         subpdt(ig,l)/imicro-pdt(ig,l)
462          ENDDO
463       ENDDO
464       
465c------ Tracers tendencies pdqcloud
466       DO l=1,nlay
467         DO ig=1,ngrid
468         
469            pdqcloudco2(ig,l,igcm_co2_ice) =
470     &        subpdq(ig,l,igcm_co2_ice)/imicro
471     &       - pdq(ig,l,igcm_co2_ice)
472            pdqcloudco2(ig,l,igcm_co2) =
473     &        subpdq(ig,l,igcm_co2)/imicro
474     &       - pdq(ig,l,igcm_co2)
475         ENDDO
476       ENDDO
477
478       
479!      call WRITEdiagfi(ngrid,"co2cloud00","co2 traceur","kg/kg",1,
480!     &      pq(1,:,igcm_co2_ice) + ptimestep*
481!     &      (pdq(1,:,igcm_co2_ice) + pdqcloudco2(1,:,igcm_co2_ice)))
482     
483       
484       IF(microphysco2) THEN
485        DO l=1,nlay
486         DO ig=1,ngrid
487            pdqcloudco2(ig,l,igcm_ccnco2_mass) =
488     &        subpdq(ig,l,igcm_ccnco2_mass)/imicro
489     &       - pdq(ig,l,igcm_ccnco2_mass)
490            pdqcloudco2(ig,l,igcm_ccnco2_number) =
491     &        subpdq(ig,l,igcm_ccnco2_number)/imicro
492     &       - pdq(ig,l,igcm_ccnco2_number)
493         ENDDO
494        ENDDO
495       ENDIF
496
497       
498       IF(scavenging) THEN
499        DO l=1,nlay
500         DO ig=1,ngrid
501            pdqcloudco2(ig,l,igcm_dust_mass) =
502     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
503     &       - pdq(ig,l,igcm_dust_mass)
504            pdqcloudco2(ig,l,igcm_dust_number) =
505     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
506     &       - pdq(ig,l,igcm_dust_number)
507         ENDDO
508        ENDDO
509        ENDIF
510
511c       ENDIF
512c------- Due to stepped entry, other processes tendencies can add up to negative values
513c------- Therefore, enforce positive values and conserve mass
514
515
516       IF(microphysco2) THEN
517        DO l=1,nlay
518         DO ig=1,ngrid
519          IF ((pq(ig,l,igcm_ccnco2_number) +
520     &      ptimestep* (pdq(ig,l,igcm_ccnco2_number) +
521     &        pdqcloudco2(ig,l,igcm_ccnco2_number))
522     &           .lt. 0)
523     &   .or. (pq(ig,l,igcm_ccnco2_mass) +
524     &      ptimestep* (pdq(ig,l,igcm_ccnco2_mass) +
525     &        pdqcloudco2(ig,l,igcm_ccnco2_mass))
526     &           .lt. 0)) THEN
527
528         pdqcloudco2(ig,l,igcm_ccnco2_number) =
529     &     - pq(ig,l,igcm_ccnco2_number)/ptimestep
530     &     - pdq(ig,l,igcm_ccnco2_number) +0
531
532         pdqcloudco2(ig,l,igcm_dust_number) = 
533     &     -pdqcloudco2(ig,l,igcm_ccnco2_number)
534
535         pdqcloudco2(ig,l,igcm_ccnco2_mass) =
536     &     - pq(ig,l,igcm_ccnco2_mass)/ptimestep
537     &     - pdq(ig,l,igcm_ccnco2_mass)+0
538
539         pdqcloudco2(ig,l,igcm_dust_mass) =
540     &     -pdqcloudco2(ig,l,igcm_ccnco2_mass)
541
542          ENDIF
543         ENDDO
544        ENDDO
545       ENDIF
546
547     
548       IF(scavenging) THEN
549          DO l=1,nlay
550             DO ig=1,ngrid
551                IF ( (pq(ig,l,igcm_dust_number) +
552     &               ptimestep* (pdq(ig,l,igcm_dust_number) +
553     &               pdqcloudco2(ig,l,igcm_dust_number)) .lt. 0.)
554     &               .or. (pq(ig,l,igcm_dust_mass)+
555     &               ptimestep* (pdq(ig,l,igcm_dust_mass) +
556     &               pdqcloudco2(ig,l,igcm_dust_mass))
557     &              .lt. 0.)) then
558
559                   pdqcloudco2(ig,l,igcm_dust_number) =
560     &                  - pq(ig,l,igcm_dust_number)/ptimestep
561     &                  - pdq(ig,l,igcm_dust_number)+0
562
563                   pdqcloudco2(ig,l,igcm_ccnco2_number) = 
564     &                  -pdqcloudco2(ig,l,igcm_dust_number)
565
566                   pdqcloudco2(ig,l,igcm_dust_mass) =
567     &                  - pq(ig,l,igcm_dust_mass)/ptimestep
568     &                  - pdq(ig,l,igcm_dust_mass) +0
569
570                   pdqcloudco2(ig,l,igcm_ccnco2_mass) =
571     &                  -pdqcloudco2(ig,l,igcm_dust_mass)
572                ENDIF
573             ENDDO
574          ENDDO
575       ENDIF !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq
576
577     
578        DO l=1,nlay
579         DO ig=1,ngrid
580          IF (pq(ig,l,igcm_co2_ice) + ptimestep*
581     &       (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice))
582     &       .lt. 1.e-25) THEN
583           pdqcloudco2(ig,l,igcm_co2_ice) =
584     &     - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)
585     &            +1.e-25
586           pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
587          ENDIF
588         ENDDO
589        ENDDO
590       
591
592
593
594c------Update the ice and dust particle size "riceco2" for output or photochemistry
595c------Only rsedcloudco2 is used for the co2 (cloud) cycle
596
597       IF(scavenging) THEN
598          DO l=1, nlay
599             DO ig=1,ngrid
600
601c        call updaterdust(
602c     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
603c     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
604c     &    pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
605c     &    pq(ig,l,igcm_dust_number) +                 ! dust number
606c     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
607c     &    pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
608c     &    rdust(ig,l))
609c         write(*,*) "in co2clouds, rdust(ig,l)= ",rdust(ig,l)
610                mdustJA= pq(ig,l,igcm_dust_mass) +             
611     &               (pdq(ig,l,igcm_dust_mass) +             
612     &               pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep
613                ndustJA=pq(ig,l,igcm_dust_number) +
614     &               (pdq(ig,l,igcm_dust_number) +
615     &               pdqcloudco2(ig,l,igcm_dust_number))*ptimestep
616               if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt.
617     &               1.e-30 *tauscaling(ig))) then
618                   rdust(ig,l)=1.e-10
619                else
620                   rdust(ig,l)=(3./4./pi/2500.*mdustJA/ndustJA)**(1./3.)
621                   rdust(ig,l)=min(rdust(ig,l),5.e-4)
622                   rdust(ig,l)=max(rdust(ig,l),1.e-10)
623                endif
624             ENDDO
625          ENDDO
626       ENDIF
627       
628       
629      IF(microphysco2) THEN
630       
631       DO l=1, nlay
632         DO ig=1,ngrid
633
634c     call updaterice_microco2(
635c     &    pq(ig,l,igcm_co2_ice) +                    ! ice mass
636c     &   (pdq(ig,l,igcm_co2_ice) +                   ! ice mass
637c     &    pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep,    ! ice mass
638c     &    pq(ig,l,igcm_ccnco2_mass) +                   ! ccn mass
639c     &   (pdq(ig,l,igcm_ccnco2_mass) +                  ! ccn mass
640c     &    pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep,   ! ccn mass
641c     &    pq(ig,l,igcm_ccnco2_number) +                 ! ccn number
642c     &   (pdq(ig,l,igcm_ccnco2_number) +                ! ccn number
643c     &    pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep, ! ccn number
644c     &    tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
645c        write(*,*) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l)
646       
647         
648        Niceco2=pq(ig,l,igcm_co2_ice) +                   
649     &       (pdq(ig,l,igcm_co2_ice) +
650     &       pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep
651        Nccnco2=max((pq(ig,l,igcm_ccnco2_number) +                 
652     &       (pdq(ig,l,igcm_ccnco2_number) +               
653     &       pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)*
654     &       tauscaling(ig),1.e-30)
655        Qccnco2=max((pq(ig,l,igcm_ccnco2_mass) +                 
656     &       (pdq(ig,l,igcm_ccnco2_mass) +               
657     &       pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)*
658     &       tauscaling(ig),1.e-30)
659        rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 + Qccnco2*rho_dust)
660     &       / (Niceco2 + Qccnco2)
661c        rhocloudco2(ig,l) = min(max(rhocloudco2t,rho_ice_co2),rho_dust)
662
663c        write(*,*) "test, nccnco2 =",nccnco22
664
665
666        riceco2(ig,l)= Niceco2*3.0/
667     &       (4.0*rho_ice_co2*pi*Nccnco2)
668     &        +rdust(ig,l)*rdust(ig,l)*rdust(ig,l)
669
670        riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
671        write(*,*) "In co2cloud, after loop, riceco2 =",riceco2(ig,l)
672        write(*,*) "In co2cloud, after loop, rhoco2 ="
673     &       ,rhocloudco2t(ig,l)
674
675        call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
676     &             tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
677
678        write(*,*) "In co2cloud, after loop and update, riceco2 ="
679     &       ,riceco2(ig,l)
680        write(*,*) "In co2cloud, after loop and update, rhoco2 ="
681     &       ,rhocloudco2t(ig,l)
682
683        if ( Niceco2
684     &       .le. 1.e-23 .or. riceco2(ig,l) .le. 1.e-10 .or.
685     &       riceco2(ig,l) .ge. 4.999e-4) then ! .or. riceco2(ig,l) .gt. 1.e-4  ) then
686           riceco2(ig,l)=0.
687           
688           !NO CLOUD : RESET TRACER AND CONSERVE MASS
689           pdqcloudco2(ig,l,igcm_co2)= pq(ig,l,igcm_co2_ice)
690     &          /ptimestep+pdq(ig,l,igcm_co2_ice)
691           
692           pdqcloudco2(ig,l,igcm_co2_ice)=-pq(ig,l,igcm_co2_ice)
693     &          /ptimestep-pdq(ig,l,igcm_co2_ice)
694           
695           pdqcloudco2(ig,l,igcm_ccnco2_mass)=
696     &          -pq(ig,l,igcm_ccnco2_mass)
697     &          /ptimestep-pdq(ig,l,igcm_ccnco2_mass)
698
699           pdqcloudco2(ig,l,igcm_ccnco2_number)=
700     &          -pq(ig,l,igcm_ccnco2_number)
701     &          /ptimestep-pdq(ig,l,igcm_ccnco2_number)
702
703           pdqcloudco2(ig,l,igcm_dust_number)=
704     &          pq(ig,l,igcm_ccnco2_number)
705     &          /ptimestep+pdq(ig,l,igcm_ccnco2_number)
706           
707           pdqcloudco2(ig,l,igcm_dust_mass)=
708     &          pq(ig,l,igcm_ccnco2_mass)
709     &          /ptimestep+pdq(ig,l,igcm_ccnco2_mass)
710c$$$
711
712c$$$           
713       endif
714     
715c       write(*,*) "in co2clouds, riceco2(ig,l)v2= ",riceco2(ig,l)
716           
717      ENDDO
718      ENDDO
719     
720      ELSE ! no microphys ! not of concern for co2 clouds  - listo
721       
722       ENDIF ! of IF(microphysco2)
723     
724     
725c     TO CHEK for relevancy - listo
726
727c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
728c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
729c     Then that should not affect the ice particle radius
730      do ig=1,ngrid
731        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
732          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
733     &    riceco2(ig,2)=riceco2(ig,3)
734          riceco2(ig,1)=riceco2(ig,2)
735        end if
736      end do
737       
738       
739       DO l=1,nlay
740         DO ig=1,ngrid
741           rsedcloudco2(ig,l)=max(riceco2(ig,l)*
742     &           (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed),
743     &                    rdust(ig,l))
744          rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-4)
745         ENDDO
746       ENDDO
747       
748       call co2sat(ngrid*nlay,pt,pplay,zqsatco2)
749         do ig=1,ngrid
750            do l=1,nlay
751               satuco2(ig,l) = pq(ig,l,igcm_co2)*
752     &              (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
753                 
754               write(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
755            enddo
756         enddo
757       call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",1,
758     &        satuco2)
759       call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
760     &        ,1,riceco2)
761! or output in diagfi.nc (for testphys1d)
762c         call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
763c         call WRITEDIAGFI(ngrid,'temp','Temperature ',
764c     &                       'K JA',1,pt)
765         
766      call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",1,
767     &   rsedcloudco2)
768     
769! used for rad. transfer calculations
770! nuice is constant because a lognormal distribution is prescribed
771c      nuice(1:ngrid,1:nlay)=nuice_ref
772
773
774
775c=======================================================================
776
777      END
778 
Note: See TracBrowser for help on using the repository browser.