source: trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F @ 1771

Last change on this file since 1771 was 1711, checked in by mvals, 8 years ago

Mars GCM

  • implementation of a sub-grid water cloud fraction scheme (by A. Pottier): If CLFvarying is set to true in callphys.def, then the sub-grid cloud fraction routine is applied in watercloud.F and aeropacity.F.
  • accordingly modified files are: watercloud.F, aeropacity.F, callradite.F, conf_phys.F, phyetat0.F, phyredem.F90, physiq_mod.F, testphys1d.F, callkeys.h, newtsart.F, start2archive.F, lect_start_archive.F
  • added file: tcondwater.F90, used by watercloud.F to calculate the condensation temperature of water
  • watercloud.F, aeropacity.F, callradite.F are converted to module files as watercloud_mod.F, aeropacity_mod.F, callradite_mod.F

MV

File size: 27.3 KB
RevLine 
[1711]1      MODULE aeropacity_mod
2
3      IMPLICIT NONE
4
5      CONTAINS
6
[38]7      SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
[1353]8     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
[1711]9     &    nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
10     &    clearsky,totcloudfrac)
[38]11                                                   
12! to use  'getin'
[1036]13      USE ioipsl_getincom, only: getin
14      use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
[1224]15     &                      igcm_dust_submicron, rho_dust, rho_ice,
16     &                      nqdust
[1543]17      use geometry_mod, only: latitude ! grid point latitudes (rad)
[1541]18      use comgeomfi_h, only: sinlat ! sines of grid point latitudes
[1375]19#ifdef DUSTSTORM
[1543]20      use geometry_mod, only: longitude
[1375]21      use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
22#endif
[1226]23      use planete_h
24      USE comcstfi_h
[1246]25      use dimradmars_mod, only: naerkind, name_iaer,
26     &            iaerdust,tauvis,
27     &            iaer_dust_conrath,iaer_dust_doubleq,
28     &            iaer_dust_submicron,iaer_h2o_ice
[38]29       IMPLICIT NONE
30c=======================================================================
31c   subject:
32c   --------
33c   Computing aerosol optical depth in each gridbox.
34c
35c   author: F.Forget
36c   ------
37c   update F. Montmessin (water ice scheme)
38c      and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
39c   update J.-B. Madeleine 2008-2009:
40c       - added 3D scattering by aerosols;
41c       - dustopacity transferred from physiq.F to callradite.F,
42c           and renamed into aeropacity.F;
[607]43c   update E. Millour, march 2012:
44c         - reference pressure is now set to 610Pa (not 700Pa)
[38]45c   
46c   input:
47c   -----
48c   ngrid             Number of gridpoint of horizontal grid
49c   nlayer            Number of layer
50c   nq                Number of tracer
51c   zday                  Date (time since Ls=0, in martian days)
52c   ls                Solar longitude (Ls) , radian
53c   pplay,pplev       pressure (Pa) in the middle and boundary of each layer
54c   pq                Dust mixing ratio (used if tracer =T and active=T).
55c   reffrad(ngrid,nlayer,naerkind)  Aerosol effective radius
[1047]56c   QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
57c   QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
58c   omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
59c   omegaREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
[38]60c
61c   output:
62c   -------
[607]63c   tauref       Prescribed mean column optical depth at 610 Pa
[38]64c   tau          Column total visible dust optical depth at each point
65c   aerosol      aerosol(ig,l,1) is the dust optical
66c                depth in layer l, grid point ig
67
68c
69c=======================================================================
70#include "callkeys.h"
71
72c-----------------------------------------------------------------------
73c
74c    Declarations :
75c    --------------
76c
77c    Input/Output
78c    ------------
79      INTEGER ngrid,nlayer,nq
80
81      REAL ls,zday,expfactor   
82      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
83      REAL pq(ngrid,nlayer,nq)
84      REAL tauref(ngrid), tau(ngrid,naerkind)
85      REAL aerosol(ngrid,nlayer,naerkind)
[1047]86      REAL dsodust(ngrid,nlayer)
[38]87      REAL reffrad(ngrid,nlayer,naerkind)
88      REAL nueffrad(ngrid,nlayer,naerkind)
[1047]89      REAL QREFvis3d(ngrid,nlayer,naerkind)
90      REAL QREFir3d(ngrid,nlayer,naerkind)
91      REAL omegaREFvis3d(ngrid,nlayer,naerkind)
92      REAL omegaREFir3d(ngrid,nlayer,naerkind)
[1711]93      real,intent(in) :: totcloudfrac(ngrid) ! total cloud fraction
94      logical,intent(in) :: clearsky ! true for part without clouds,
95                                     ! false for part with clouds (total or sub-grid clouds)
96      real :: CLFtot ! total cloud fraction
[38]97c
98c    Local variables :
99c    -----------------
100      INTEGER l,ig,iq,i,j
101      INTEGER iaer           ! Aerosol index
[1047]102      real topdust(ngrid)
[38]103      real zlsconst, zp
104      real taueq,tauS,tauN
105c     Mean Qext(vis)/Qext(ir) profile
[1047]106      real msolsir(nlayer,naerkind)
[38]107c     Mean Qext(ir)/Qabs(ir) profile
[1047]108      real mqextsqabs(nlayer,naerkind)
[38]109c     Variables used when multiple particle sizes are used
110c       for dust or water ice particles in the radiative transfer
111c       (see callradite.F for more information).
[1047]112      REAL taudusttmp(ngrid)! Temporary dust opacity
[38]113                               !   used before scaling
[1047]114      REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
115      REAL taudustvis(ngrid) ! Dust opacity after scaling
116      REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
[38]117                               !   "seen" by the GCM.
[1047]118      REAL taucloudvis(ngrid)! Cloud opacity at visible
[38]119                               !   reference wavelength
[1047]120      REAL taucloudtes(ngrid)! Cloud opacity at infrared
[38]121                               !   reference wavelength using
122                               !   Qabs instead of Qext
123                               !   (direct comparison with TES)
[1224]124      REAL topdust0(ngrid)
[83]125
[1375]126#ifdef DUSTSTORM
127!! Local dust storms
128      logical localstorm        ! =true to create a local dust storm
129      real taulocref,ztoploc,radloc,lonloc,latloc  ! local dust storm parameters
130      real reffstorm, yeah
131      REAL ray(ngrid) ! distance from dust storm center
132      REAL tauuser(ngrid) ! opacity perturbation due to dust storm
133      REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
134      REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
135      real l_top ! layer of the storm's top
136      REAL zalt(ngrid, nlayer) ! useful factor to compute l_top
137#endif
138
[38]139c   local saved variables
140c   ---------------------
141
142c     Level under which the dust mixing ratio is held constant
143c       when computing the dust opacity in each layer
144c       (this applies when doubleq and active are true)
[1208]145      INTEGER, PARAMETER :: cstdustlevel0 = 7
146      INTEGER, SAVE      :: cstdustlevel
[38]147
[607]148      LOGICAL,SAVE :: firstcall=.true.
[38]149
150! indexes of water ice and dust tracers:
151      INTEGER,SAVE :: i_ice=0  ! water ice
[607]152      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
[38]153      CHARACTER(LEN=20) :: txt ! to temporarly store text
154      CHARACTER(LEN=1) :: txt2 ! to temporarly store text
155! indexes of dust scatterers:
156      INTEGER,SAVE :: naerdust ! number of dust scatterers
157
158      tau(1:ngrid,1:naerkind)=0
159
160! identify tracers
161
162      IF (firstcall) THEN
163        ! identify scatterers that are dust
164        naerdust=0
165        DO iaer=1,naerkind
166          txt=name_iaer(iaer)
167          IF (txt(1:4).eq."dust") THEN
168            naerdust=naerdust+1
169            iaerdust(naerdust)=iaer
170          ENDIF
171        ENDDO
172        ! identify tracers which are dust
173        i=0
174        DO iq=1,nq
175          txt=noms(iq)
176          IF (txt(1:4).eq."dust") THEN
177          i=i+1
178          nqdust(i)=iq
179          ENDIF
180        ENDDO
181
182        IF (water.AND.activice) THEN
183          i_ice=igcm_h2o_ice
184          write(*,*) "aeropacity: i_ice=",i_ice
185        ENDIF
186
187c       typical profile of solsir and (1-w)^(-1):
188        msolsir(1:nlayer,1:naerkind)=0
189        mqextsqabs(1:nlayer,1:naerkind)=0
[222]190        WRITE(*,*) "Typical profiles of Qext(vis)/Qext(IR)"
191        WRITE(*,*) "  and Qext(IR)/Qabs(IR):"
[38]192        DO iaer = 1, naerkind ! Loop on aerosol kind
193          WRITE(*,*) "Aerosol # ",iaer
194          DO l=1,nlayer
[1047]195            DO ig=1,ngrid
[38]196              msolsir(l,iaer)=msolsir(l,iaer)+
197     &              QREFvis3d(ig,l,iaer)/
198     &              QREFir3d(ig,l,iaer)
199              mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
200     &              (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
201            ENDDO
[1047]202            msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
203            mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
[38]204          ENDDO
205          WRITE(*,*) "solsir: ",msolsir(:,iaer)
206          WRITE(*,*) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
207        ENDDO
208
209!       load value of tauvis from callphys.def (if given there,
210!       otherwise default value read from starfi.nc file will be used)
211        call getin("tauvis",tauvis)
212
[1208]213        IF (freedust) THEN
214          cstdustlevel = 1
215        ELSE
216          cstdustlevel = cstdustlevel0
217        ENDIF
218
219
[1375]220#ifndef DUSTSTORM
[38]221        firstcall=.false.
[1375]222#endif
[38]223
224      END IF
225
[607]226c     Vertical column optical depth at "odpref" Pa
227c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1088]228      IF(freedust) THEN
229         tauref(:) = 0. ! tauref is computed after, instead of being forced
230
231      ELSE IF(iaervar.eq.1) THEN
[1047]232         do ig=1, ngrid
[38]233          tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
234                                       ! or read in starfi
235        end do
236      ELSE IF (iaervar.eq.2) THEN   ! << "Viking" Scenario>>
237
238        tauref(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1
239        do ig=2,ngrid
240          tauref(ig) = tauref(1)
241        end do
242
243      ELSE IF (iaervar.eq.3) THEN  ! << "MGS" scenario >>
244
245        taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls-4.363)))**14
246        tauS= 0.1 +(0.5-0.1)  *(cos(0.5*(ls-4.363)))**14
247        tauN = 0.1
248c          if (peri_day.eq.150) then
249c            tauS=0.1
250c            tauN=0.1 +(0.5-0.1)  *(cos(0.5*(ls+pi-4.363)))**14
251c            taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls+pi-4.363)))**14
252c           endif
[1047]253        do ig=1,ngrid
[1541]254          if (latitude(ig).ge.0) then
[1047]255          ! Northern hemisphere
256            tauref(ig)= tauN +
[1541]257     &      (taueq-tauN)*0.5*(1+tanh((45-latitude(ig)*180./pi)*6/60))
[1047]258          else
259          ! Southern hemisphere
260            tauref(ig)= tauS +
[1541]261     &      (taueq-tauS)*0.5*(1+tanh((45+latitude(ig)*180./pi)*6/60))
[1047]262          endif
263        enddo ! of do ig=1,ngrid
[38]264      ELSE IF (iaervar.eq.5) THEN   ! << Escalier Scenario>>
265c         tauref(1) = 0.2
266c         if ((ls.ge.210.*pi/180.).and.(ls.le.330.*pi/180.))
267c    &                              tauref(1) = 2.5
268        tauref(1) = 2.5
269        if ((ls.ge.30.*pi/180.).and.(ls.le.150.*pi/180.))
270     &                              tauref(1) = .2
271
272        do ig=2,ngrid
273          tauref(ig) = tauref(1)
274        end do
[1278]275      ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
276      ! clim, cold or warm synthetic scenarios
[677]277        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
[1502]278      ELSE IF ((iaervar.ge.24).and.(iaervar.le.32))
[607]279     &     THEN  ! << MY... dust scenarios >>
280        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
281      ELSE IF ((iaervar.eq.4).or.
282     &         ((iaervar.ge.124).and.(iaervar.le.126))) THEN
283       ! "old" TES assimation dust scenario (values at 700Pa in files!)
284        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
[38]285      ELSE
286        stop 'problem with iaervar in aeropacity.F'
287      ENDIF
288
289c -----------------------------------------------------------------
290c Computing the opacity in each layer
291c -----------------------------------------------------------------
292
293      DO iaer = 1, naerkind ! Loop on aerosol kind
294c     --------------------------------------------
295        aerkind: SELECT CASE (name_iaer(iaer))
296c==================================================================
297        CASE("dust_conrath") aerkind      ! Typical dust profile
298c==================================================================
299
300c       Altitude of the top of the dust layer
301c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
302        zlsconst=SIN(ls-2.76)
303        if (iddist.eq.1) then
304          do ig=1,ngrid
305             topdust(ig)=topdustref         ! constant dust layer top
306          end do
307
308        else if (iddist.eq.2) then          ! "Viking" scenario
309          do ig=1,ngrid
[1224]310            ! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
311            ! in the Viking year scenario
312            topdust0(ig)=60. -22.*sinlat(ig)**2
[38]313            topdust(ig)=topdust0(ig)+18.*zlsconst
314          end do
315
316        else if(iddist.eq.3) then         !"MGS" scenario
317          do ig=1,ngrid
318            topdust(ig)=60.+18.*zlsconst
[1541]319     &                -(32+18*zlsconst)*sin(latitude(ig))**4
320     &                 - 8*zlsconst*(sin(latitude(ig)))**5
[38]321          end do
322        endif
323
324c       Optical depth in each layer :
325c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
326        if(iddist.ge.1) then
327
328          expfactor=0.
329          DO l=1,nlayer
330            DO ig=1,ngrid
331c             Typical mixing ratio profile
[607]332              if(pplay(ig,l).gt.odpref
[38]333     $                        /(988.**(topdust(ig)/70.))) then
[607]334                zp=(odpref/pplay(ig,l))**(70./topdust(ig))
[38]335                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
336              else   
337                expfactor=1.e-3
338              endif
339c             Vertical scaling function
340              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
341     &          expfactor *
342     &          QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
343            ENDDO
344          ENDDO
345
346        else if(iddist.eq.0) then   
347c         old dust vertical distribution function (pollack90)
348          DO l=1,nlayer
349             DO ig=1,ngrid
[607]350                zp=odpref/pplay(ig,l)
351                aerosol(ig,l,1)= tauref(ig)/odpref *
[38]352     s           (pplev(ig,l)-pplev(ig,l+1))
353     s           *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 )
354             ENDDO
355          ENDDO
356        end if
357
358c==================================================================
359        CASE("dust_doubleq") aerkind! Two-moment scheme for dust
360c        (transport of mass and number mixing ratio)
361c==================================================================
362
363          DO l=1,nlayer
364            IF (l.LE.cstdustlevel) THEN
365c           Opacity in the first levels is held constant to
366c             avoid unrealistic values due to constant lifting:
367              DO ig=1,ngrid
368                aerosol(ig,l,iaer) =
369     &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
370     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
371     &          pq(ig,cstdustlevel,igcm_dust_mass) *
372     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
[1353]373                ! DENSITY SCALED OPACITY IN INFRARED:
374                dsodust(ig,l) =
375     &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
376     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
377     &          pq(ig,cstdustlevel,igcm_dust_mass)
[38]378              ENDDO
379            ELSE
380              DO ig=1,ngrid
381                aerosol(ig,l,iaer) =
382     &          (  0.75 * QREFvis3d(ig,l,iaer) /
383     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
384     &          pq(ig,l,igcm_dust_mass) *
385     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
[1353]386                ! DENSITY SCALED OPACITY IN INFRARED:
387                dsodust(ig,l) =
388     &          (  0.75 * QREFir3d(ig,l,iaer) /
389     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
390     &          pq(ig,l,igcm_dust_mass)
[38]391              ENDDO
392            ENDIF
393          ENDDO
394
395c==================================================================
396        CASE("dust_submicron") aerkind   ! Small dust population
397c==================================================================
398
399          DO l=1,nlayer
400            IF (l.LE.cstdustlevel) THEN
401c           Opacity in the first levels is held constant to
402c             avoid unrealistic values due to constant lifting:
403              DO ig=1,ngrid
404                aerosol(ig,l,iaer) =
405     &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
406     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
407     &          pq(ig,cstdustlevel,igcm_dust_submicron) *
408     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
409              ENDDO
410            ELSE
411              DO ig=1,ngrid
412                aerosol(ig,l,iaer) =
413     &          (  0.75 * QREFvis3d(ig,l,iaer) /
414     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
415     &          pq(ig,l,igcm_dust_submicron) *
416     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
417              ENDDO
418            ENDIF
419          ENDDO
420
421c==================================================================
422        CASE("h2o_ice") aerkind             ! Water ice crystals
423c==================================================================
424
425c       1. Initialization
426        aerosol(1:ngrid,1:nlayer,iaer) = 0.
427        taucloudvis(1:ngrid) = 0.
428        taucloudtes(1:ngrid) = 0.
429c       2. Opacity calculation
[1711]430        ! NO CLOUDS
431        IF (clearsky) THEN
432            aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
433        ! CLOUDSs
434        ELSE ! else (clearsky)
435          DO ig=1, ngrid
436            DO l=1,nlayer
437              aerosol(ig,l,iaer) = max(1E-20,
438     &          (  0.75 * QREFvis3d(ig,l,iaer) /
439     &          ( rho_ice * reffrad(ig,l,iaer) )  ) *
440     &          pq(ig,l,i_ice) *
441     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
[38]442     &                              )
[1711]443              taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
444              taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
445     &          QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
446     &          ( 1.E0 - omegaREFir3d(ig,l,iaer) )
447            ENDDO
[38]448          ENDDO
[1711]449          ! SUB-GRID SCALE CLOUDS
450          IF (CLFvarying) THEN
451             DO ig=1, ngrid
452                DO l=1,nlayer-1
453                   CLFtot  = max(totcloudfrac(ig),0.01)
454                   aerosol(ig,l,iaer)=
455     &                    aerosol(ig,l,iaer)/CLFtot
456                   aerosol(ig,l,iaer) =
457     &                    max(aerosol(ig,l,iaer),1.e-9)
458                ENDDO
459             ENDDO
460!          ELSE ! else (CLFvarying)
461!             DO ig=1, ngrid
462!                DO l=1,nlayer-1 ! to stop the rad tran bug
463!                   CLFtot  = CLFfixval
464!                   aerosol(ig,l,iaer)=
465!     &                    aerosol(ig,l,iaer)/CLFtot
466!                   aerosol(ig,l,iaer) =
467!     &                    max(aerosol(ig,l,iaer),1.e-9)
468!                ENDDO
469!             ENDDO
470          ENDIF ! end (CLFvarying)             
471        ENDIF ! end (clearsky)
472
[520]473c       3. Outputs -- Now done in physiq.F
474!        IF (ngrid.NE.1) THEN
[1047]475!          CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
[520]476!     &      ' ',2,taucloudvis)
[1047]477!          CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
[520]478!     &      ' ',2,taucloudtes)
479!          IF (callstats) THEN
[1047]480!            CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
[520]481!     &        ' ',2,taucloudvis)
[1047]482!            CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
[520]483!     &        ' ',2,taucloudtes)
484!          ENDIF
485!        ELSE
486!         CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
487!        ENDIF
[38]488c==================================================================
489        END SELECT aerkind
490c     -----------------------------------
491      ENDDO ! iaer (loop on aerosol kind)
492
493c -----------------------------------------------------------------
494c Rescaling each layer to reproduce the choosen (or assimilated)
495c   dust extinction opacity at visible reference wavelength, which
[607]496c   is originally scaled to an equivalent odpref Pa pressure surface.
[38]497c -----------------------------------------------------------------
498
[1375]499#ifdef DUSTSTORM
500c -----------------------------------------------------------------
[1410]501! Calculate reference opacity without perturbation
[1375]502c -----------------------------------------------------------------
503      IF (firstcall) THEN
504        DO iaer=1,naerdust
505          DO l=1,nlayer
506            DO ig=1,ngrid
507              tauref(ig) = tauref(ig) +
508     &                    aerosol(ig,l,iaerdust(iaer))
509            ENDDO
510          ENDDO
511        ENDDO
512        tauref(:) = tauref(:) * odpref / pplev(:,1)
[1410]513
[1375]514c--------------------------------------------------
[1410]515c Get parameters of the opacity perturbation
[1375]516c--------------------------------------------------
[1410]517        iaer=1  ! just change dust
[1375]518
519        write(*,*) "Add a local storm ?"
520        localstorm=.true. ! default value
521        call getin("localstorm",localstorm)
522        write(*,*) " localstorm = ",localstorm
523
524        IF (localstorm) THEN
525          WRITE(*,*) "********************"
526          WRITE(*,*) "ADDING A LOCAL STORM"
527          WRITE(*,*) "********************"
528
529          write(*,*) "ref opacity of local dust storm"
530              taulocref = 4.25 ! default value
531              call getin("taulocref",taulocref)
532              write(*,*) " taulocref = ",taulocref
533
534          write(*,*) "target altitude of local storm (km)"
535              ztoploc = 10.0 ! default value
536              call getin("ztoploc",ztoploc)
537              write(*,*) " ztoploc = ",ztoploc
538
539          write(*,*) "radius of dust storm (degree)"
540              radloc = 0.5 ! default value
541              call getin("radloc",radloc)
542              write(*,*) " radloc = ",radloc
543
544          write(*,*) "center longitude of storm (deg)"
545              lonloc = 25.0 ! default value
546              call getin("lonloc",lonloc)
547              write(*,*) " lonloc = ",lonloc
548
549          write(*,*) "center latitude of storm (deg)"
550              latloc = -2.5 ! default value
551              call getin("latloc",latloc)
552              write(*,*) " latloc = ",latloc
553       
554          write(*,*) "reff storm (mic) 0. for background"
555              reffstorm = 0.0 ! default value
556              call getin("reffstorm",reffstorm)
557              write(*,*) " reffstorm = ",reffstorm
558
[1410]559!! LOOP: modify opacity
[1375]560      DO ig=1,ngrid
561
[1410]562      !! distance to the center:
[1541]563      ray(ig)=SQRT((latitude(ig)*180./pi-latloc)**2 +
564     &          (longitude(ig)*180./pi -lonloc)**2)
[1375]565
566      !! transition factor for storm
[1410]567      !! factor is hardcoded -- increase it to steepen
[1375]568      yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
569
[1410]570      !! new opacity field
571      !! -- add an opacity set to taulocref
572      !! -- the additional reference opacity will
573      !!      thus be taulocref*odpref/pplev
574      tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref ,
575     &          taulocref * yeah )
[1375]576
[1410]577      !! compute l_top
[1375]578          DO l=1,nlayer
579            zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
580     &                      / g / 44.01
581     &                    * 8.31 * 210.
582                IF (     (ztoploc .lt. zalt(ig,l)  )
583     &          .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
584          ENDDO
585
[1410]586     !! change reffrad if ever needed
[1375]587      IF (reffstorm .gt. 0.) THEN
588          DO l=1,nlayer
589             IF (l .lt. l_top+1) THEN
590                reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
591     &          * 1.e-6 * yeah )
592             ENDIF
593          ENDDO
594      ENDIF
595
[1410]596      ENDDO
597!! END LOOP
[1375]598
[1410]599      !! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
[1375]600      DO ig=1,ngrid
601          int_factor(ig)=0.
602          DO l=1,nlayer
603             IF (l .lt. l_top+1) THEN
604                      int_factor(ig) =
605     &                int_factor(ig) +
606     &          (  0.75 * QREFvis3d(ig,l,iaer) /
607     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
608     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
609             ENDIF
610          ENDDO
611          DO l=1, nlayer
[1410]612          !! Mass mixing ratio perturbation due to local dust storm in each layer
[1375]613          more_dust(ig,l,1)=
614     &                     (tauuser(ig)-(tauref(ig)
615     &                      * pplev(ig,1) /odpref)) /
616     &                      int_factor(ig)
617          more_dust(ig,l,2)=
618     &                     (tauuser(ig)-(tauref(ig) *
619     &                     pplev(ig,1) /odpref))
620     &                      / int_factor(ig) *
621     &                     ((ref_r0/reffrad(ig,l,iaer))**3)
622     &                      * r3n_q
623          ENDDO
624      ENDDO
625
[1410]626      !! quantity of dust for each layer with the addition of the perturbation
627      DO l=1, l_top
[1376]628          pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
[1410]629     .            + more_dust(:,l,1)
[1376]630          pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
[1410]631     .            + more_dust(:,l,2)
632      ENDDO
633      ENDIF !! IF (localstorm)
634      tauref(:)=0.
635      ENDIF !! IF (firstcall)
[1375]636#endif
637
[1088]638      IF (freedust) THEN
639          tauscaling(:) = 1.
640
641      ELSE
642c       Temporary scaling factor
643        taudusttmp(1:ngrid)=0.
644        DO iaer=1,naerdust
645          DO l=1,nlayer
646            DO ig=1,ngrid
647c             Scaling factor
648              taudusttmp(ig) = taudusttmp(ig) +
649     &                         aerosol(ig,l,iaerdust(iaer))
650            ENDDO
[38]651          ENDDO
652        ENDDO
[358]653
[1088]654c       Saved scaling factor
655        DO ig=1,ngrid
656            tauscaling(ig) = tauref(ig) *
657     &                       pplev(ig,1) / odpref / taudusttmp(ig)
658        ENDDO
[358]659
[1410]660      ENDIF ! IF (freedust)
[1088]661
[358]662c     Opacity computation
[38]663      DO iaer=1,naerdust
664        DO l=1,nlayer
665          DO ig=1,ngrid
666            aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
[411]667     &                aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
[38]668          ENDDO
669        ENDDO
670      ENDDO
[1088]671
[1353]672      DO l=1,nlayer
673        DO ig=1,ngrid
674          dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
675        ENDDO
676      ENDDO
677
[1088]678      IF (freedust) THEN
[1208]679        ! tauref has been initialized to 0 before.
680        DO iaer=1,naerdust
681          DO l=1,nlayer
682            DO ig=1,ngrid
[1410]683#ifdef DUSTSTORM
684      !! recalculate opacity because storm perturbation has been added
685      IF (firstcall) THEN
686              aerosol(ig,l,iaer) =
687     &          (  0.75 * QREFvis3d(ig,l,iaer) /
688     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
689     &          pq(ig,l,igcm_dust_mass) *
690     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
691      ENDIF
692#endif
[1208]693              tauref(ig) = tauref(ig) +
694     &                    aerosol(ig,l,iaerdust(iaer))
695            ENDDO
[1088]696          ENDDO
697        ENDDO
[1208]698        tauref(:) = tauref(:) * odpref / pplev(:,1)
[1088]699      ENDIF
[411]700     
701c output for debug
[420]702c        IF (ngrid.NE.1) THEN
[1047]703c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
[420]704c     &      '#',2,taudusttmp)
[1047]705c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
[420]706c     &      '#',2,tauscaling)
707c        ELSE
[1047]708c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
[420]709c     &      '#',0,taudusttmp)
[1047]710c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
[420]711c     &      '#',0,tauscaling)
712c        ENDIF
[38]713c -----------------------------------------------------------------
714c Column integrated visible optical depth in each point
715c -----------------------------------------------------------------
716      DO iaer=1,naerkind
717        do l=1,nlayer
718           do ig=1,ngrid
719             tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
720           end do
721        end do
722      ENDDO
[1375]723
724#ifdef DUSTSTORM
725      IF (firstcall) THEN
726        firstcall=.false.
727      ENDIF
728#endif
729
[38]730c -----------------------------------------------------------------
731c Density scaled opacity and column opacity output
732c -----------------------------------------------------------------
[420]733c      dsodust(1:ngrid,1:nlayer) = 0.
734c      DO iaer=1,naerdust
[1047]735c        DO l=1,nlayer
[420]736c          DO ig=1,ngrid
737c            dsodust(ig,l) = dsodust(ig,l) +
738c     &                      aerosol(ig,l,iaerdust(iaer)) * g /
739c     &                      (pplev(ig,l) - pplev(ig,l+1))
740c          ENDDO
741c        ENDDO
742c        IF (ngrid.NE.1) THEN
743c          write(txt2,'(i1.1)') iaer
[1047]744c          call WRITEDIAGFI(ngrid,'taudust'//txt2,
[420]745c     &                    'Dust col opacity',
746c     &                    ' ',2,tau(1,iaerdust(iaer)))
747c          IF (callstats) THEN
[1047]748c            CALL wstats(ngrid,'taudust'//txt2,
[420]749c     &                 'Dust col opacity',
750c     &                 ' ',2,tau(1,iaerdust(iaer)))
751c          ENDIF
752c        ENDIF
753c      ENDDO
[38]754
[420]755c      IF (ngrid.NE.1) THEN
[1047]756c       CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
[38]757c    &                    'm2.kg-1',3,dsodust)
[420]758c        IF (callstats) THEN
[1047]759c          CALL wstats(ngrid,'dsodust',
[420]760c     &               'tau*g/dp',
761c     &               'm2.kg-1',3,dsodust)
762c        ENDIF
763c      ELSE
764c        CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
765c     &                   dsodust)
766c      ENDIF ! of IF (ngrid.NE.1)
[38]767c -----------------------------------------------------------------
[1711]768
769      END SUBROUTINE aeropacity
770     
771      END MODULE aeropacity_mod
Note: See TracBrowser for help on using the repository browser.