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

Last change on this file since 1775 was 1754, checked in by aslmd, 8 years ago

LMDZ.MARS : tcondco2 is not present, prevents the model to compile. commented the line.
  • Property svn:executable set to *
File size: 44.7 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,
6     &                rsedcloud,rhocloud,pzlev,pdqs_sedco2)
7! to use  'getin'
8      use dimradmars_mod, only: naerkind
9      USE comcstfi_h
10      USE ioipsl_getincom
11      USE updaterad
12      use conc_mod, only: mmean
13      use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
14     &     igcm_dust_mass, igcm_dust_number,igcm_h2o_ice,
15     &     igcm_ccn_mass,igcm_ccn_number,
16     &     igcm_ccnco2_mass, igcm_ccnco2_number,
17     &     rho_dust, nuiceco2_sed, nuiceco2_ref,
18     &     rho_ice_co2,r3n_q,rho_ice,nuice_sed
19     
20
21      IMPLICIT NONE
22
23
24#include "datafile.h"
25#include "callkeys.h"
26#include "microphys.h"
27
28
29c=======================================================================
30c CO2 clouds formation
31c
32c  There is a time loop specific to cloud formation
33c  due to timescales smaller than the GCM integration timestep.
34c  microphysics subroutine is improvedCO2clouds.F
35c 
36c  The co2 clouds tracers (co2_ice, ccn mass and concentration) are
37c  sedimented at each microtimestep. pdqs_sedco2 keeps track of the
38c  CO2 flux at the surface
39c
40c  Authors: 09/2016 Joachim Audouard & Constantino Listowski
41c  Adaptation of the water ice clouds scheme (with specific microphysics)
42c  of Montmessin, Navarro & al.
43c
44c 07/2017 J.Audouard
45c Several logicals can be set to .true. in callphys.def
46c co2clouds=.true. call this routine
47c co2useh2o=.true. allow the use of water ice particles as CCN for CO2
48c meteo_flux=.true. supply meteoritic particles
49c CLFvaryingCO2=.true. allows a subgrid temperature distribution
50c of amplitude spantCO2(=integer in callphys.def)
51c
52c=======================================================================
53
54c-----------------------------------------------------------------------
55c   declarations:
56c   -------------
57
58c   Inputs:
59c   ------
60
61      INTEGER ngrid,nlay
62      INTEGER nq                 ! nombre de traceurs
63      REAL ptimestep            ! pas de temps physique (s)
64      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
65      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
66      REAL pdpsrf(ngrid)         ! tendence surf pressure
67      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
68      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
69      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
70      real,intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
71
72      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
73      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
74
75      real rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
76                                ! used for nucleation of CO2 on ice-coated ccns
77      DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
78      REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
79      REAL tauscaling(ngrid)   ! Convertion factor for dust amount
80      real rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
81
82c   Outputs:
83c   -------
84
85      real pdqcloudco2(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
86      REAL pdtcloudco2(ngrid,nlay)    ! tendence temperature due
87                                   ! a la chaleur latente
88
89      DOUBLE PRECISION riceco2(ngrid,nlay)    ! Ice mass mean radius (m)
90                               ! (r_c in montmessin_2004)
91      REAL nuice(ngrid,nlay)   ! Estimated effective variance
92                               !   of the size distribution
93      real rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
94      real rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
95      real rhocloudco2t(ngrid,nlay)  ! Cloud density (kg.m-3)
96      real  pdqs_sedco2(ngrid) ! CO2 flux at the surface
97c   local:
98c   ------
99      !water
100      real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
101      real rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
102      ! for ice radius computation
103     
104      REAl ccntyp
105     
106      ! for time loop
107      INTEGER microstep  ! time subsampling step variable
108      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
109      SAVE imicro
110      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
111      SAVE microtimestep
112     
113      ! tendency given by clouds (inside the micro loop)
114      REAL subpdqcloudco2(ngrid,nlay,nq) ! cf. pdqcloud
115      REAL subpdtcloudco2(ngrid,nlay)    ! cf. pdtcloud
116
117      ! global tendency (clouds+physics)
118      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
119      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
120      real wq(ngrid,nlay+1)  !  ! displaced tracer mass (kg.m-2) during microtimestep because sedim (?/m-2)
121
122      REAL satuco2(ngrid,nlay)  ! co2 satu ratio for output
123      REAL zqsatco2(ngrid,nlay) ! saturation co2
124
125      INTEGER iq,ig,l,i
126      LOGICAL,SAVE :: firstcall=.true.
127      DOUBLE PRECISION Nccnco2, Niceco2,Nco2,mdustJA,ndustJA
128      DOUBLE PRECISION Qccnco2
129      real :: beta
130
131      real epaisseur (ngrid,nlay) ! Layer thickness (m)
132      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
133      double precision diff,diff0
134      real tempo_traceur_t(ngrid,nlay)
135      real tempo_traceurs(ngrid,nlay,nq)
136      real sav_trac(ngrid,nlay,nq)
137      real pdqsed(ngrid,nlay,nq)
138      REAL lw                   !Latent heat of sublimation (J.kg-1)
139      REAL,save :: l0,l1,l2,l3,l4 !Coeffs for lw
140
141      DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) ! Nb particules intégré
142      DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:)
143      DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:)
144      DOUBLE PRECISION :: myT
145!     What we need for Qext reading and tau computation
146      DOUBLE PRECISION vrat_cld ! Volume ratio
147      DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
148      SAVE rb_cldco2
149      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-11 ! Minimum radius (m)
150      DOUBLE PRECISION, PARAMETER :: rmax_cld = 3.e-6 ! Maximum radius (m)
151      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-12
152! Minimum boundary radius (m)
153      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m)
154
155      DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
156      DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
157      REAL sigma_iceco2 ! Variance of the ice and CCN distributions
158
159      logical :: file_ok
160      double precision :: radv(10000),Qextv1mic(10000)
161      double precision :: Qext1bins(100),Qtemp
162      double precision :: ltemp1(10000),ltemp2(10000)
163      integer :: nelem,lebon1,lebon2,uQext
164      save Qext1bins
165      character(len=100) scanline 
166      DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2
167      DOUBLE PRECISION Qext1bins2(ngrid,nlay)   
168      DOUBLE PRECISION tau1mic(ngrid) !co2 ice column opacity at 1µm
169     
170        ! For sub grid T distribution
171
172      REAL zt(ngrid,nlay)       ! local value of temperature
173      REAL :: zq(ngrid, nlay,nq)
174
175      REAL :: tcond(ngrid,nlay)
176      REAL ::  zqvap(ngrid,nlay)
177      REAL :: zqice(ngrid,nlay)
178      REAL ::  spant,zdelt ! delta T for the temperature distribution
179      REAL :: zteff(ngrid, nlay)! effective temperature in the cloud,neb
180      REAL :: pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud
181      REAL :: cloudfrac(ngrid,nlay) ! cloud fraction
182      REAL :: mincloud ! min cloud frac
183c      logical :: CLFvaryingCO2
184c     ** un petit test de coherence
185c       --------------------------
186
187      IF (firstcall) THEN
188        if (nq.gt.nqmx) then
189           write(*,*) 'stop in co2cloud (nq.gt.nqmx)!'
190           write(*,*) 'nq=',nq,' nqmx=',nqmx
191           stop
192        endif
193        write(*,*) "co2cloud.F: rho_ice_co2 = ",rho_ice_co2
194        write(*,*) "co2cloud: igcm_co2=",igcm_co2
195        write(*,*) "            igcm_co2_ice=",igcm_co2_ice
196               
197        write(*,*) "time subsampling for microphysic ?"
198#ifdef MESOSCALE
199        imicro = 2
200#else
201        imicro = 30
202#endif
203        call getin("imicro",imicro)
204        write(*,*)"imicro = ",imicro
205       
206        microtimestep = ptimestep/real(imicro)
207        write(*,*)"Physical timestep is",ptimestep
208        write(*,*)"CO2 Microphysics timestep is",microtimestep
209     
210        ! Values for latent heat computation
211        l0=595594d0     
212        l1=903.111d0   
213        l2=-11.5959d0   
214        l3=0.0528288d0
215        l4=-0.000103183d0
216c$$$       
217c$$$        if (meteo_flux_number .ne. 0) then
218c$$$
219c$$$        write(*,*) "Warning ! Meteoritic flux of dust is turned on"
220c$$$        write(*,*) "Dust mass = ",meteo_flux_mass
221c$$$        write(*,*) "Dust number = ",meteo_flux_number
222c$$$        write(*,*) "Are added at the z-level (km)",meteo_alt
223c$$$        write(*,*) "Every timestep in co2cloud.F"
224c$$$         endif
225         
226        if (.not. allocated(memdMMccn)) allocate(memdMMccn(ngrid,nlay))
227        if (.not. allocated(memdNNccn)) allocate(memdNNccn(ngrid,nlay))
228        if (.not. allocated(memdMMh2o)) allocate(memdMMh2o(ngrid,nlay))
229       
230        memdMMccn(:,:)=0.
231        memdMMh2o(:,:)=0.
232        memdNNccn(:,:)=0.
233
234c     Compute the size bins of the distribution of CO2 ice particles
235c --> used for opacity calculations
236
237c       rad_cldco2 is the primary radius grid used for microphysics computation.
238c       The grid spacing is computed assuming a constant volume ratio
239c       between two consecutive bins; i.e. vrat_cld.
240c       vrat_cld is determined from the boundary values of the size grid:
241c       rmin_cld and rmax_cld.
242c       The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
243c       dr_cld is the width of each rad_cldco2 bin.
244        sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
245c       Volume ratio between two adjacent bins
246   !     vrat_cld
247        vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
248        vrat_cld = exp(vrat_cld)
249        rb_cldco2(1)  = rbmin_cld
250        rad_cldco2(1) = rmin_cld
251        vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
252        do i=1,nbinco2_cld-1
253          rad_cldco2(i+1)  = rad_cldco2(i) * vrat_cld**(1./3.)
254          vol_cld(i+1)  = vol_cld(i) * vrat_cld
255        enddo
256        do i=1,nbinco2_cld
257          rb_cldco2(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
258     &      rad_cldco2(i)
259          dr_cld(i)  = rb_cldco2(i+1) - rb_cldco2(i)
260        enddo
261        rb_cldco2(nbinco2_cld+1) = rbmax_cld
262        dr_cld(nbinco2_cld)   = rb_cldco2(nbinco2_cld+1) -
263     &       rb_cldco2(nbinco2_cld)
264
265c   read the Qext values
266        INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//
267     &       '/optprop_co2ice_1mic.dat', EXIST=file_ok)
268        IF (.not. file_ok) THEN
269           write(*,*) 'file optprop_co2ice_1mic.dat should be in '
270     &          ,datafile
271           STOP
272        endif
273        open(newunit=uQext,file=trim(datafile)//
274     &       '/optprop_co2ice_1mic.dat'
275     &       ,FORM='formatted')
276        read(uQext,*) !skip 1 line
277        do i=1,10000
278           read(uQext,'(E11.5)') radv(i)
279        enddo
280        read(uQext,*) !skip 1 line
281        do i=1,10000
282           read(uQext,'(E11.5)') Qextv1mic(i)
283        enddo
284        close(uQext)
285c   innterpol the Qext values
286        !rice_out=rad_cldco2
287        do i=1,nbinco2_cld
288           ltemp1=abs(radv(:)-rb_cldco2(i))
289           ltemp2=abs(radv(:)-rb_cldco2(i+1))
290           lebon1=minloc(ltemp1,DIM=1)
291           lebon2=minloc(ltemp2,DIM=1)
292           nelem=lebon2-lebon1+1.
293           Qtemp=0d0
294           do l=0,nelem
295              Qtemp=Qtemp+Qextv1mic(lebon1+l) !mean value in the interval
296           enddo
297           Qtemp=Qtemp/nelem
298           Qext1bins(i)=Qtemp
299        enddo
300        Qext1bins(:)=Qext1bins(:)*rad_cldco2(:)*rad_cldco2(:)*pi
301!     The actuall tau computation and output is performed in co2cloud.F
302
303        print*,'--------------------------------------------'
304        print*,'Microphysics co2: size bin-Qext information:'
305        print*,'   i, rad_cldco2(i), Qext1bins(i)'
306        do i=1,nbinco2_cld
307          write(*,'(i3,3x,3(e12.6,4x))') i, rad_cldco2(i),
308     &    Qext1bins(i)
309        enddo
310        print*,'--------------------------------------------'
311
312
313        do i=1,nbinco2_cld+1
314            rb_cldco2(i) = log(rb_cldco2(i))  !! we save that so that it is not computed  at each timestep and gridpoint
315        enddo
316        if (CLFvaryingCO2) then
317          write(*,*)
318          write(*,*) "CLFvaryingCO2 is set to true is callphys.def"
319          write(*,*) "The temperature field is enlarged to +/-",spantCO2
320          write(*,*) "for the CO2 microphysics "
321        endif
322        firstcall=.false.
323      ENDIF                     ! of IF (firstcall)
324   
325c-----Initialization
326      dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
327      beta=0.85
328      subpdq(1:ngrid,1:nlay,1:nq) = 0
329      subpdt(1:ngrid,1:nlay)      = 0
330      subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
331      subpdtcloudco2(1:ngrid,1:nlay)      = 0
332c$$$
333c$$$       call WRITEDIAGFI(ngrid,"pzlay","pzlay","km",3,
334c$$$     &        pzlay)
335c$$$       call WRITEDIAGFI(ngrid,"pplay","pplay","Pa",3,
336c$$$     &        pplay)
337
338      wq(:,:)=0
339      ! default value if no ice
340      rhocloudco2(1:ngrid,1:nlay) = rho_dust
341      rhocloudco2t(1:ngrid,1:nlay) = rho_dust
342      epaisseur(1:ngrid,1:nlay)=0
343      masse(1:ngrid,1:nlay)=0
344
345      sav_trac(1:ngrid,1:nlay,1:nq)=0
346      pdqsed(1:ngrid,1:nlay,1:nq)=0
347     
348      do  l=1,nlay
349        do ig=1, ngrid
350          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
351          epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
352       enddo
353      enddo
354
355      !CLFvaryingCO2=.true.
356     
357c-------------------------------------------------------------------
358c   0.  Representation of sub-grid water ice clouds
359c------------------
360      IF (CLFvaryingCO2) THEN
361         spant=spantCO2         ! delta T for the temprature distribution
362         mincloud=0.1           ! min cloudfrac when there is ice 
363        ! IF (flagcloudco2) THEN
364        !    write(*,*) "CLFCO2 Delta T", spant
365        !    write(*,*) "CLFCO2 mincloud", mincloud
366        !    flagcloudco2=.false.
367        ! END IF
368         
369c-----Tendencies
370         DO l=1,nlay
371            DO ig=1,ngrid
372               zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
373            ENDDO
374         ENDDO
375         DO l=1,nlay
376            DO ig=1,ngrid
377               DO iq=1,nq
378                  zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
379               ENDDO
380            ENDDO
381         ENDDO
382         zqvap=zq(:,:,igcm_co2)
383         zqice=zq(:,:,igcm_co2_ice)
384!! AS : this routine is not present?
385!         CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
386! A tester:            CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
387         zdelt=spant 
388         DO l=1,nlay
389            DO ig=1,ngrid
390               IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt)) THEN !Toute la fraction est saturée
391                  zteff(ig,l)=zt(ig,l)
392                  cloudfrac(ig,l)=1.
393               ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN !Rien n'est saturé
394                  zteff(ig,l)=zt(ig,l)-zdelt
395                  cloudfrac(ig,l)=mincloud
396               ELSE
397                  cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/
398     &                 (2.0*zdelt)
399                  zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Ja Temperature moyenne de la fraction nuageuse
400                 
401               END IF
402               zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep
403               IF (cloudfrac(ig,l).le. mincloud) THEN
404                  cloudfrac(ig,l)=mincloud
405               ELSE IF (cloudfrac(ig,l).ge. 1) THEN
406                  cloudfrac(ig,l)=1.
407               END IF
408            ENDDO
409         ENDDO
410!     Totalcloud frac of the column missing here
411c-----------------------
412c-----No sub-grid cloud representation (CLFvarying=false)
413      ELSE
414         DO l=1,nlay
415            DO ig=1,ngrid
416               zteff(ig,l)=pt(ig,l)
417            END DO
418         END DO
419      END IF                    ! end if (CLFvaryingco2)--------------------------------------------
420c   1.  Tendencies:
421c------------------
422 
423c------ Effective tracers quantities in the cloud fraction
424        IF (CLFvaryingCO2) THEN     
425           pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
426   
427c           pqeff(:,:,igcm_ccn_mass) =pq(:,:,igcm_ccn_mass)/
428c     &                                cloudfrac(:,:)
429c           pqeff(:,:,igcm_ccn_number)=
430c     &     pq(:,:,igcm_ccn_number)/cloudfrac(:,:)
431c           pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)/
432c     &                               cloudfrac(:,:)
433           pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/
434     &                                cloudfrac(:,:)
435           pqeff(:,:,igcm_ccnco2_number)=
436     &     pq(:,:,igcm_ccnco2_number)/cloudfrac(:,:)
437           pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/
438     &                               cloudfrac(:,:)
439
440        ELSE
441           pqeff(:,:,:)=pq(:,:,:)
442c           pqeff(:,:,igcm_ccn_mass)= pq(:,:,igcm_ccn_mass)
443c           pqeff(:,:,igcm_ccn_number)= pq(:,:,igcm_ccn_number)
444c           pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)
445c           pqeff(:,:,igcm_ccnco2_mass)= pq(:,:,igcm_ccnco2_mass)
446c           pqeff(:,:,igcm_ccnco2_number)= pq(:,:,igcm_ccnco2_number)
447c           pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)
448        END IF     
449        tempo_traceur_t(1:ngrid,1:nlay)=zteff(1:ngrid,1:nlay)
450      tempo_traceurs(1:ngrid,1:nlay,1:nq)=pqeff(1:ngrid,1:nlay,1:nq)
451
452c------------------------------------------------------------------
453c Time subsampling for microphysics
454c------------------------------------------------------------------
455      DO microstep=1,imicro
456c------ Temperature tendency subpdt
457        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
458        ! If imicro=1 subpdt is the same as pdt
459        DO l=1,nlay
460          DO ig=1,ngrid
461
462             subpdt(ig,l) = subpdt(ig,l)
463     &            + pdt(ig,l)   ! At each micro timestep we add pdt in order to have a stepped entry
464             subpdq(ig,l,igcm_dust_mass) =
465     &            subpdq(ig,l,igcm_dust_mass)
466     &            + pdq(ig,l,igcm_dust_mass)
467             subpdq(ig,l,igcm_dust_number) =
468     &            subpdq(ig,l,igcm_dust_number)
469     &            + pdq(ig,l,igcm_dust_number)
470             subpdq(ig,l,igcm_ccnco2_mass) =
471     &            subpdq(ig,l,igcm_ccnco2_mass)
472     &            + pdq(ig,l,igcm_ccnco2_mass)
473             subpdq(ig,l,igcm_ccnco2_number) =
474     &            subpdq(ig,l,igcm_ccnco2_number)
475     &            + pdq(ig,l,igcm_ccnco2_number)
476             subpdq(ig,l,igcm_co2_ice) =
477     &            subpdq(ig,l,igcm_co2_ice)
478     &            + pdq(ig,l,igcm_co2_ice)
479             subpdq(ig,l,igcm_co2) =
480     &            subpdq(ig,l,igcm_co2)
481     &            + pdq(ig,l,igcm_co2)
482             subpdq(ig,l,igcm_h2o_ice) =
483     &            subpdq(ig,l,igcm_h2o_ice)
484     &            + pdq(ig,l,igcm_h2o_ice)
485             subpdq(ig,l,igcm_ccn_mass) =
486     &            subpdq(ig,l,igcm_ccn_mass)
487     &            + pdq(ig,l,igcm_ccn_mass)
488             subpdq(ig,l,igcm_ccn_number) =
489     &            subpdq(ig,l,igcm_ccn_number)
490     &            + pdq(ig,l,igcm_ccn_number)
491          ENDDO
492       ENDDO
493
494c add meteoritic flux of dust (old version)
495cNew version (John Plane values) are added in improvedCO2clouds
496    !convert meteo_alt (in km) to z-level
497        !pzlay altitudes of the layers
498c$$$!zlev altitudes (nlayl+1) of the boundaries
499c$$$       if (meteo_flux_number .ge. 0) then
500c$$$          do ig=1, ngrid
501c$$$             l=1
502c$$$             do  while ( (((meteo_alt .ge. pplev(ig,l)) .and.
503c$$$     &            (meteo_alt .le. pplev(ig,l+1))) .eq. .false.)
504c$$$     &            .and. (l .lt. nlay) )
505c$$$                l=l+1
506c$$$             enddo
507c$$$             meteo_lvl=l
508c$$$             subpdq(ig,meteo_lvl,igcm_dust_mass)=
509c$$$     &            subpdq(ig,meteo_lvl,igcm_dust_mass)
510c$$$     &            +meteo_flux_mass
511c$$$             subpdq(ig,meteo_lvl,igcm_dust_number)=
512c$$$     &            subpdq(ig,meteo_lvl,igcm_dust_number)
513c$$$     &            +meteo_flux_number
514c$$$          enddo
515c$$$       endif
516c-------------------------------------------------------------------
517c   2.  Main call to the cloud schemes:
518c------------------------------------------------
519      CALL improvedCO2clouds(ngrid,nlay,microtimestep,
520     &     pplay,pzlev,zteff,subpdt,
521     &     pqeff,subpdq,subpdqcloudco2,subpdtcloudco2,
522     &     nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)
523c-------------------------------------------------------------------
524c   3.  Updating tendencies after cloud scheme:
525c-----------------------------------------------
526          DO l=1,nlay
527            DO ig=1,ngrid
528               subpdt(ig,l) =
529     &              subpdt(ig,l) + subpdtcloudco2(ig,l)
530               subpdq(ig,l,igcm_dust_mass) =
531     &              subpdq(ig,l,igcm_dust_mass)
532     &              + subpdqcloudco2(ig,l,igcm_dust_mass)
533               subpdq(ig,l,igcm_dust_number) =
534     &              subpdq(ig,l,igcm_dust_number)
535     &              + subpdqcloudco2(ig,l,igcm_dust_number)
536               subpdq(ig,l,igcm_ccnco2_mass) =
537     &              subpdq(ig,l,igcm_ccnco2_mass)
538     &              + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
539               subpdq(ig,l,igcm_ccnco2_number) =
540     &              subpdq(ig,l,igcm_ccnco2_number)
541     &              + subpdqcloudco2(ig,l,igcm_ccnco2_number)
542               subpdq(ig,l,igcm_co2_ice) =
543     &              subpdq(ig,l,igcm_co2_ice)
544     &              + subpdqcloudco2(ig,l,igcm_co2_ice)
545               subpdq(ig,l,igcm_co2) =
546     &              subpdq(ig,l,igcm_co2)
547     &              + subpdqcloudco2(ig,l,igcm_co2)
548               subpdq(ig,l,igcm_h2o_ice) =
549     &              subpdq(ig,l,igcm_h2o_ice)
550     &              + subpdqcloudco2(ig,l,igcm_h2o_ice)
551               subpdq(ig,l,igcm_ccn_mass) =
552     &              subpdq(ig,l,igcm_ccn_mass)
553     &              + subpdqcloudco2(ig,l,igcm_ccn_mass)
554               subpdq(ig,l,igcm_ccn_number) =
555     &              subpdq(ig,l,igcm_ccn_number)
556     &              + subpdqcloudco2(ig,l,igcm_ccn_number)
557            ENDDO
558         ENDDO
559       
560!Sedimentation
561!Update traceurs, compute rsed
562       
563       DO l=1, nlay
564          DO ig=1,ngrid             
565             tempo_traceur_t(ig,l)=zteff(ig,l)+subpdt(ig,l)
566     &            *microtimestep
567             tempo_traceurs(ig,l,:)=pqeff(ig,l,:)
568     &            +subpdq(ig,l,:)*microtimestep
569
570             rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*
571     &            tempo_traceur_t(ig,l)-2.87e-6*
572     &            tempo_traceur_t(ig,l)*tempo_traceur_t(ig,l))
573             
574             rho_ice_co2=rho_ice_co2T(ig,l)
575             Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30)
576             Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number),
577     &            1.e-30)
578             Qccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_mass),
579     &            1.e-30)
580             mdustJA= tempo_traceurs(ig,l,igcm_dust_mass)
581             ndustJA=tempo_traceurs(ig,l,igcm_dust_number)
582             if ((Nccnco2 .lt. tauscaling(ig)) .or. (Qccnco2 .lt.
583     &            1.e-30 *tauscaling(ig))) then
584                rdust(ig,l)=1.e-10
585             else
586                rdust(ig,l)=(3./4./pi/2500.*Qccnco2/Nccnco2)**(1./3.)
587                rdust(ig,l)=max(rdust(ig,l),1.e-10)
588c     rdust(ig,l)=min(rdust(ig,l),5.e-6)   
589             end if
590             rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2
591     &            + Qccnco2*tauscaling(ig)*rho_dust)
592     &            / (Niceco2 + Qccnco2*tauscaling(ig))
593             
594             rhocloudco2t(ig,l)=min(max(rhocloudco2t(ig,l)
595     &            ,rho_ice_co2),rho_dust)
596             riceco2(ig,l)=(Niceco2*3.0/
597     &            (4.0*rho_ice_co2*pi*Nccnco2
598     &            *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
599     &            *rdust(ig,l))**(1.0/3.0)
600             riceco2(ig,l)=max(1.e-10,riceco2(ig,l))
601             
602             rsedcloudco2(ig,l)=max(riceco2(ig,l)*
603     &            (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed),
604     &            rdust(ig,l))
605             
606          ENDDO
607       ENDDO
608       
609!     Gravitational sedimentation
610       
611       sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
612       sav_trac(:,:,igcm_ccnco2_mass)=
613     &      tempo_traceurs(:,:,igcm_ccnco2_mass)
614       sav_trac(:,:,igcm_ccnco2_number)=
615     &      tempo_traceurs(:,:,igcm_ccnco2_number)
616       
617      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
618     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
619     &     rsedcloudco2,rhocloudco2t,
620     &     tempo_traceurs(:,:,igcm_co2_ice),wq,beta) !  3 traceurs
621     
622!     sedim at the surface of co2 ice : keep track of it for physiq_mod
623      do ig=1,ngrid
624         pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep
625      end do
626
627      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
628     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
629     &     rsedcloudco2,rhocloudco2t,
630     &     tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta)
631     
632      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
633     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
634     &     rsedcloudco2,rhocloudco2t,
635     &     tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta)
636           
637      DO l = 1, nlay !Compute tendencies
638         DO ig=1,ngrid
639            pdqsed(ig,l,igcm_ccnco2_mass)=
640     &           (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
641     &           sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
642            pdqsed(ig,l,igcm_ccnco2_number)=
643     &           (tempo_traceurs(ig,l,igcm_ccnco2_number)-
644     &           sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
645            pdqsed(ig,l,igcm_co2_ice)=
646     &           (tempo_traceurs(ig,l,igcm_co2_ice)-
647     &           sav_trac(ig,l,igcm_co2_ice))/microtimestep
648         ENDDO
649      ENDDO
650            !pdqsed est la tendance due a la sedimentation
651     
652      DO l=1,nlay
653         DO ig=1,ngrid
654            subpdq(ig,l,igcm_ccnco2_mass) =
655     &           subpdq(ig,l,igcm_ccnco2_mass)
656     &           +pdqsed(ig,l,igcm_ccnco2_mass)
657            subpdq(ig,l,igcm_ccnco2_number) =
658     &           subpdq(ig,l,igcm_ccnco2_number)
659     &           +pdqsed(ig,l,igcm_ccnco2_number)
660            subpdq(ig,l,igcm_co2_ice) =
661     &           subpdq(ig,l,igcm_co2_ice)
662     &           +pdqsed(ig,l,igcm_co2_ice)
663         ENDDO
664      ENDDO   
665 
666      ENDDO                     ! of DO microstep=1,imicro
667     
668c-------------------------------------------------------------------
669c   6.  Compute final tendencies after time loop:
670c------------------------------------------------
671c CO2 flux at surface (kg.m-2.s-1)
672      do ig=1,ngrid
673         pdqs_sedco2(ig)=pdqs_sedco2(ig)/real(imicro)
674      enddo
675
676c------ Temperature tendency pdtcloud
677       DO l=1,nlay
678         DO ig=1,ngrid
679             pdtcloudco2(ig,l) =
680     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
681          ENDDO
682       ENDDO
683
684c------ Tracers tendencies pdqcloud
685       DO l=1,nlay
686         DO ig=1,ngrid
687         
688            pdqcloudco2(ig,l,igcm_co2_ice) =
689     &        subpdq(ig,l,igcm_co2_ice)/real(imicro)
690     &       - pdq(ig,l,igcm_co2_ice)
691            pdqcloudco2(ig,l,igcm_co2) =
692     &        subpdq(ig,l,igcm_co2)/real(imicro)
693     &       - pdq(ig,l,igcm_co2)
694            pdqcloudco2(ig,l,igcm_h2o_ice) =
695     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro)
696     &       - pdq(ig,l,igcm_h2o_ice)
697         ENDDO
698       ENDDO
699
700        DO l=1,nlay
701         DO ig=1,ngrid
702            pdqcloudco2(ig,l,igcm_ccnco2_mass) =
703     &        subpdq(ig,l,igcm_ccnco2_mass)/real(imicro)
704     &       - pdq(ig,l,igcm_ccnco2_mass)
705            pdqcloudco2(ig,l,igcm_ccnco2_number) =
706     &        subpdq(ig,l,igcm_ccnco2_number)/real(imicro)
707     &       - pdq(ig,l,igcm_ccnco2_number)
708            pdqcloudco2(ig,l,igcm_ccn_mass) =
709     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro)
710     &       - pdq(ig,l,igcm_ccn_mass)
711            pdqcloudco2(ig,l,igcm_ccn_number) =
712     &        subpdq(ig,l,igcm_ccn_number)/real(imicro)
713     &       - pdq(ig,l,igcm_ccn_number)
714         ENDDO
715        ENDDO
716
717       
718        DO l=1,nlay
719         DO ig=1,ngrid
720            pdqcloudco2(ig,l,igcm_dust_mass) =
721     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
722     &       - pdq(ig,l,igcm_dust_mass)
723            pdqcloudco2(ig,l,igcm_dust_number) =
724     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
725     &       - pdq(ig,l,igcm_dust_number)
726         ENDDO
727        ENDDO
728
729c------- Due to stepped entry, other processes tendencies can add up to negative values
730c------- Therefore, enforce positive values and conserve mass
731
732
733
734        DO l=1,nlay
735           DO ig=1,ngrid
736              IF ((pqeff(ig,l,igcm_ccnco2_number) +
737     &             ptimestep* (pdq(ig,l,igcm_ccnco2_number) +
738     &             pdqcloudco2(ig,l,igcm_ccnco2_number))
739     &             .lt. 1.e-20)
740     &             .or. (pqeff(ig,l,igcm_ccnco2_mass) +
741     &             ptimestep* (pdq(ig,l,igcm_ccnco2_mass) +
742     &             pdqcloudco2(ig,l,igcm_ccnco2_mass))
743     &             .lt. 1.e-30)) THEN
744                 
745                 pdqcloudco2(ig,l,igcm_ccnco2_number) =
746     &                - pqeff(ig,l,igcm_ccnco2_number)/ptimestep
747     &                - pdq(ig,l,igcm_ccnco2_number)+1.e-20
748                 pdqcloudco2(ig,l,igcm_dust_number) = 
749     &                -pdqcloudco2(ig,l,igcm_ccnco2_number)
750                 pdqcloudco2(ig,l,igcm_ccnco2_mass) =
751     &                - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep
752     &                - pdq(ig,l,igcm_ccnco2_mass)+1.e-30
753                 pdqcloudco2(ig,l,igcm_dust_mass) =
754     &                -pdqcloudco2(ig,l,igcm_ccnco2_mass)
755
756              ENDIF
757           ENDDO
758        ENDDO
759       
760
761     
762        DO l=1,nlay
763           DO ig=1,ngrid
764              IF ( (pqeff(ig,l,igcm_dust_number) +
765     &             ptimestep* (pdq(ig,l,igcm_dust_number) +
766     &             pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-30)
767     &             .or. (pqeff(ig,l,igcm_dust_mass)+
768     &             ptimestep* (pdq(ig,l,igcm_dust_mass) +
769     &             pdqcloudco2(ig,l,igcm_dust_mass))
770     &             .le. 1.e-30)) then
771                 
772                 pdqcloudco2(ig,l,igcm_dust_number) =
773     &                - pqeff(ig,l,igcm_dust_number)/ptimestep
774     &                - pdq(ig,l,igcm_dust_number)+1.e-30
775                 pdqcloudco2(ig,l,igcm_ccnco2_number) = 
776     &                  -pdqcloudco2(ig,l,igcm_dust_number)
777                 pdqcloudco2(ig,l,igcm_dust_mass) =
778     &                  - pqeff(ig,l,igcm_dust_mass)/ptimestep
779     &                  - pdq(ig,l,igcm_dust_mass) +1.e-30
780                 pdqcloudco2(ig,l,igcm_ccnco2_mass) =
781     &                -pdqcloudco2(ig,l,igcm_dust_mass)
782
783              ENDIF
784           ENDDO
785        ENDDO
786        !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq
787c$$$
788c$$$     
789c$$$        DO l=1,nlay
790c$$$         DO ig=1,ngrid
791c$$$          IF (pq(ig,l,igcm_co2_ice) + ptimestep*
792c$$$     &       (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice))
793c$$$     &       .lt. 1.e-30) THEN
794c$$$           pdqcloudco2(ig,l,igcm_co2_ice) =
795c$$$     &     - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)
796c$$$           pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
797c$$$           !write(*,*) "WARNING CO2 ICE in co2cloud.F"
798c$$$
799c$$$          ENDIF   
800c$$$          IF (pq(ig,l,igcm_co2) + ptimestep*
801c$$$     &       (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2))
802c$$$     &       .lt. 0.5) THEN
803c$$$           pdqcloudco2(ig,l,igcm_co2) =
804c$$$     &     - pdq(ig,l,igcm_co2_ice) !- pdq(ig,l,igcm_co2)
805c$$$c           pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2)
806c$$$           pdqcloudco2(ig,l,igcm_co2_ice) =
807c$$$     &     - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)     
808c$$$          ! write(*,*) "WARNING CO2 in co2cloud.F"
809c$$$          ENDIF
810c$$$         ENDDO
811c$$$        ENDDO
812c$$$ 
813       
814        DO l=1, nlay
815           DO ig=1,ngrid
816
817c     call updaterice_microco2(
818c     &    pq(ig,l,igcm_co2_ice) +                    ! ice mass
819c     &   (pdq(ig,l,igcm_co2_ice) +                   ! ice mass
820c     &    pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep,    ! ice mass
821c     &    pq(ig,l,igcm_ccnco2_mass) +                   ! ccn mass
822c     &   (pdq(ig,l,igcm_ccnco2_mass) +                  ! ccn mass
823c     &    pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep,   ! ccn mass
824c     &    pq(ig,l,igcm_ccnco2_number) +                 ! ccn number
825c     &   (pdq(ig,l,igcm_ccnco2_number) +                ! ccn number
826c     &    pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep, ! ccn number
827c     &    tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
828c        write(*,*) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l)
829              Niceco2=pqeff(ig,l,igcm_co2_ice) +                   
830     &             (pdq(ig,l,igcm_co2_ice) +
831     &             pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep
832              Nco2=pqeff(ig,l,igcm_co2) +                   
833     &             (pdq(ig,l,igcm_co2) +
834     &             pdqcloudco2(ig,l,igcm_co2))*ptimestep
835              Nccnco2=max((pqeff(ig,l,igcm_ccnco2_number) +                 
836     &             (pdq(ig,l,igcm_ccnco2_number) +               
837     &             pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)*
838     &             tauscaling(ig),1.e-30)
839              Qccnco2=max((pqeff(ig,l,igcm_ccnco2_mass) +                 
840     &             (pdq(ig,l,igcm_ccnco2_mass) +               
841     &             pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)*
842     &             tauscaling(ig),1.e-30)
843             
844              if (Nccnco2 .lt. 0.1) then
845                 rdust(ig,l)=1.e-10
846              else
847       
848                 rdust(ig,l)= Qccnco2
849     &                *0.75/pi/rho_dust
850     &                / Nccnco2
851                 rdust(ig,l)= rdust(ig,l)**(1./3.)           
852                 rdust(ig,l)=max(1.e-10,rdust(ig,l))
853c     rdust(ig,l)=min(5.e-6,rdust(ig,l))
854              endif
855              myT=zteff(ig,l)+(pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
856              rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*
857     &             myT-2.87e-6* myT* myT)
858              rho_ice_co2=rho_ice_co2T(ig,l)
859
860              lw = l0 + l1 * myT + l2 *myT**2 +
861     &             l3 * myT**3 + l4 * myT**4 !J.kg-1
862
863              riceco2(ig,l)=(Niceco2*3.0/
864     &             (4.0*rho_ice_co2*pi*Nccnco2)
865     &             +rdust(ig,l)*rdust(ig,l)
866     &             *rdust(ig,l) )**(1.0/3.0)
867c    &      .or. (riceco2(ig,l) .le. rdust(ig,l))
868              if ( (Niceco2 .le. 1.e-25).or. (Nccnco2 .le. 0.1) )THEN !anciennement 200
869c           riceco2(ig,l)=0.
870
871c     &       .or. (Nccnco2 .le. 1.)
872c        endif
873!Flux chaleur latente <0 quand sublimation 
874
875           pdtcloudco2(ig,l)= pdtcloudco2(ig,l)-Niceco2*lw/cpp/ptimestep
876c$$$  !NO CLOUD : RESET TRACERS AND CONSERVE MASS
877c           if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+
878c     &          pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then
879c              pdqcloudco2(ig,l,igcm_co2)=0.
880c              pdqcloudco2(ig,l,igcm_co2_ice)=0.
881c           else
882             pdqcloudco2(ig,l,igcm_co2_ice)=-1.*pqeff(ig,l,igcm_co2_ice)
883     &             /ptimestep-pdq(ig,l,igcm_co2_ice)
884              pdqcloudco2(ig,l,igcm_co2)=-1.*
885     &             pdqcloudco2(ig,l,igcm_co2_ice)
886c           endif
887!     Reverse h2o ccn and ices into h2o tracers
888           if (memdMMccn(ig,l) .gt. 0) then
889              pdqcloudco2(ig,l,igcm_ccn_mass)=memdMMccn(ig,l)/ptimestep
890           else
891              memdMMccn(ig,l)=0.
892              pdqcloudco2(ig,l,igcm_ccn_mass)=0.
893           endif
894           if (memdNNccn(ig,l) .gt. 0) then
895             pdqcloudco2(ig,l,igcm_ccn_number)=memdNNccn(ig,l)/ptimestep
896           else
897              memdNNccn(ig,l)=0.
898              pdqcloudco2(ig,l,igcm_ccn_number)=0.
899           endif
900           if (memdMMh2o(ig,l) .gt. 0) then
901              pdqcloudco2(ig,l,igcm_h2o_ice)=memdMMh2o(ig,l)/ptimestep
902           else
903              memdMMh2o(ig,l)=0.
904              pdqcloudco2(ig,l,igcm_h2o_ice)=0.
905           endif
906c           if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+
907c     &          pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep
908c     &          .le. 1.e-30) then
909c              pdqcloudco2(ig,l,igcm_ccnco2_mass)=0.
910c              pdqcloudco2(ig,l,igcm_ccnco2_number)=0.
911c              pdqcloudco2(ig,l,igcm_co2)=0.
912c              pdqcloudco2(ig,l,igcm_co2_ice)=0.
913c           else
914              pdqcloudco2(ig,l,igcm_ccnco2_mass)=
915     &             -1.*pqeff(ig,l,igcm_ccnco2_mass)
916     &             /ptimestep-pdq(ig,l,igcm_ccnco2_mass)
917c           endif
918c           if (pq(ig,l,igcm_ccnco2_number)+
919c     &          (pdq(ig,l,igcm_ccnco2_number)+
920c     &          pdqcloudco2(ig,l,igcm_ccnco2_number))
921c     &          *ptimestep.le. 1.e-30)
922c     &          then
923c              pdqcloudco2(ig,l,igcm_ccnco2_mass)=0.
924c              pdqcloudco2(ig,l,igcm_ccnco2_number)=0.
925c              pdqcloudco2(ig,l,igcm_co2)=0.
926c              pdqcloudco2(ig,l,igcm_co2_ice)=0.
927c           else
928              pdqcloudco2(ig,l,igcm_ccnco2_number)=
929     &             -1.*pqeff(ig,l,igcm_ccnco2_number)
930     &             /ptimestep-pdq(ig,l,igcm_ccnco2_number)
931c           endif
932c           if (pq(ig,l,igcm_dust_number)+
933c     &          (pdq(ig,l,igcm_dust_number)+
934c     &          pdqcloudco2(ig,l,igcm_dust_number))
935c     &          *ptimestep.le. 1.e-30)
936c     &          then
937c              pdqcloudco2(ig,l,igcm_dust_number)=0.
938c              pdqcloudco2(ig,l,igcm_dust_mass)=0.
939c           else
940              pdqcloudco2(ig,l,igcm_dust_number)=
941     &             pqeff(ig,l,igcm_ccnco2_number)
942     &             /ptimestep+pdq(ig,l,igcm_ccnco2_number)
943     &             -memdNNccn(ig,l)/ptimestep
944c           endif
945c           if (pq(ig,l,igcm_dust_mass)+
946c     &          (pdq(ig,l,igcm_dust_mass)+
947c     &          pdqcloudco2(ig,l,igcm_dust_mass))
948c     &          *ptimestep .le. 1.e-30)
949c     &          then
950c              pdqcloudco2(ig,l,igcm_dust_number)=0.
951c              pdqcloudco2(ig,l,igcm_dust_mass)=0.
952c           else
953              pdqcloudco2(ig,l,igcm_dust_mass)=
954     &             pqeff(ig,l,igcm_ccnco2_mass)
955     &             /ptimestep+pdq(ig,l,igcm_ccnco2_mass)
956     &             -(memdMMh2o(ig,l)+memdMMccn(ig,l))/ptimestep
957c           endif
958           memdMMccn(ig,l)=0.
959           memdMMh2o(ig,l)=0.
960           memdNNccn(ig,l)=0.
961           riceco2(ig,l)=0.
962        endif
963c     Compute opacities
964      No=Nccnco2
965      Rn=-log(riceco2(ig,l))
966      n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
967      Qext1bins2(ig,l)=0.
968      do i = 1, nbinco2_cld !Qext below 50 is negligible
969         n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
970         n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
971         n_aer(i) = n_aer(i) + 0.5 * No * n_derf
972         Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)
973      enddo
974      !l'opacité de la case ig est la somme sur l de Qext1bins2
975   
976   
977      !update rice water
978        call updaterice_micro(
979     &    pqeff(ig,l,igcm_h2o_ice) +                    ! ice mass
980     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
981     &    pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
982     &    pqeff(ig,l,igcm_ccn_mass) +                   ! ccn mass
983     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
984     &    pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
985     &    pqeff(ig,l,igcm_ccn_number) +                 ! ccn number
986     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
987     &    pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
988     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
989         
990
991        call updaterdust(
992     &    pqeff(ig,l,igcm_dust_mass) +                   ! dust mass
993     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
994     &    pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
995     &    pqeff(ig,l,igcm_dust_number) +                 ! dust number
996     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
997     &    pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
998     &    rdust(ig,l))
999
1000         ENDDO
1001       ENDDO
1002       tau1mic(:)=0.
1003       do l=1,nlay
1004          do ig=1,ngrid
1005             tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
1006          enddo
1007      enddo
1008
1009c$$$
1010c$$$      if (riceco2(725,22) .ge. 1.e-10) then
1011c$$$       
1012c$$$         write(*,*) 'DIAGJA co2 ',pqeff(725,22,igcm_co2) +                   
1013c$$$     &        (pdq(725,22,igcm_co2) +
1014c$$$     &        pdqcloudco2(725,22,igcm_co2))*ptimestep
1015c$$$         write(*,*) 'DIAGJA co2_ice',pqeff(725,22,igcm_co2_ice) +
1016c$$$     &        (pdq(725,22,igcm_co2_ice) +
1017c$$$     &        pdqcloudco2(725,22,igcm_co2_ice))*ptimestep
1018c$$$         
1019c$$$         write(*,*) 'DIAGJA riceco2',riceco2(725,22)
1020c$$$         write(*,*) 'DIAGJA T',zteff(725,22) +
1021c$$$     &        (pdt(725,22) + pdtcloudco2(725,22))*ptimestep
1022c$$$         write(*,*) 'DIAG pdtcloud',pdtcloudco2(725,22)*ptimestep
1023c$$$         write(*,*) 'DIAGJA ccnNco2',pqeff(725,22,igcm_ccnco2_number)+
1024c$$$     &        (pdq(725,22,igcm_ccnco2_number) +
1025c$$$     &        pdqcloudco2(725,22,igcm_ccnco2_number))*ptimestep
1026c$$$         
1027c$$$         write(*,*) 'DIAGJA dustN',pqeff(725,22,igcm_dust_number) +
1028c$$$     &        (pdq(725,22,igcm_dust_number) +
1029c$$$     &        pdqcloudco2(725,22,igcm_dust_number))*ptimestep
1030c$$$      ENDIF
1031c$$$     
1032
1033
1034     
1035c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
1036c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1037c     Then that should not affect the ice particle radius
1038      do ig=1,ngrid
1039        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
1040          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
1041     &    riceco2(ig,2)=riceco2(ig,3)
1042          riceco2(ig,1)=riceco2(ig,2)
1043        end if
1044      end do
1045       
1046       
1047       DO l=1,nlay
1048         DO ig=1,ngrid
1049           rsedcloud(ig,l)=max(rice(ig,l)*
1050     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
1051     &                    rdust(ig,l))
1052!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
1053         ENDDO
1054      ENDDO
1055       
1056       DO l=1,nlay
1057         DO ig=1,ngrid
1058           rsedcloudco2(ig,l)=max(riceco2(ig,l)*
1059     &           (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed),
1060     &                    rdust(ig,l))
1061c          rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)
1062         ENDDO
1063       ENDDO
1064       
1065       call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep
1066     &      ,pplay,zqsatco2)
1067       do l=1,nlay
1068          do ig=1,ngrid
1069             satuco2(ig,l) = (pqeff(ig,l,igcm_co2)  +
1070     &            (pdq(ig,l,igcm_co2) +
1071     &            pdqcloudco2(ig,l,igcm_co2))*ptimestep)*
1072     &            (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
1073             !write(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
1074            enddo
1075         enddo
1076!Tout ce qui est modifié par la microphysique de CO2 doit être rapporté a cloudfrac
1077         IF (CLFvaryingCO2) THEN
1078            DO l=1,nlay
1079               DO ig=1,ngrid
1080                  pdqcloudco2(ig,l,igcm_ccn_mass)=
1081     &             pdqcloudco2(ig,l,igcm_ccn_mass)*cloudfrac(ig,l)
1082                  pdqcloudco2(ig,l,igcm_ccnco2_mass)=
1083     &             pdqcloudco2(ig,l,igcm_ccnco2_mass)*cloudfrac(ig,l)
1084                  pdqcloudco2(ig,l,igcm_ccn_number)=
1085     &             pdqcloudco2(ig,l,igcm_ccn_number)*cloudfrac(ig,l)
1086                  pdqcloudco2(ig,l,igcm_ccnco2_number)=
1087     &             pdqcloudco2(ig,l,igcm_ccnco2_number)*cloudfrac(ig,l)
1088                  pdqcloudco2(ig,l,igcm_dust_mass)=
1089     &             pdqcloudco2(ig,l,igcm_dust_mass)*cloudfrac(ig,l)
1090                  pdqcloudco2(ig,l,igcm_dust_number)=
1091     &             pdqcloudco2(ig,l,igcm_dust_number)*cloudfrac(ig,l)
1092                  pdqcloudco2(ig,l,igcm_h2o_ice)=
1093     &             pdqcloudco2(ig,l,igcm_h2o_ice)*cloudfrac(ig,l)
1094                  pdqcloudco2(ig,l,igcm_co2_ice)=
1095     &             pdqcloudco2(ig,l,igcm_co2_ice)*cloudfrac(ig,l)
1096                  pdqcloudco2(ig,l,igcm_co2)=
1097     &             pdqcloudco2(ig,l,igcm_co2)*cloudfrac(ig,l)
1098                  pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*cloudfrac(ig,l)
1099               ENDDO
1100            ENDDO   
1101         ENDIF
1102
1103         call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
1104     &        satuco2)
1105         call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
1106     &        ,3,riceco2)
1107!     or output in diagfi.nc (for testphys1d)
1108c         call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
1109c         call WRITEDIAGFI(ngrid,'temp','Temperature ',
1110c     &                       'K JA',1,pt)
1111         
1112      call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3,
1113     &   rsedcloudco2)
1114     
1115      call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"," ",2,
1116     &   tau1mic)
1117      call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"," ",3,
1118     &   cloudfrac)
1119! used for rad. transfer calculations
1120! nuice is constant because a lognormal distribution is prescribed
1121c      nuice(1:ngrid,1:nlay)=nuice_ref
1122
1123
1124
1125c=======================================================================
1126
1127      END
1128 
Note: See TracBrowser for help on using the repository browser.