source: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90 @ 799

Last change on this file since 799 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
File size: 13.4 KB
Line 
1      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
2         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
3
4       use radinc_h, only : L_TAUMAX,naerkind
5       use aerosol_mod
6       USE comgeomfi_h
7       USE tracer_h
8                 
9       implicit none
10
11!==================================================================
12!     
13!     Purpose
14!     -------
15!     Compute aerosol optical depth in each gridbox.
16!     
17!     Authors
18!     -------
19!     F. Forget
20!     F. Montmessin (water ice scheme)
21!     update J.-B. Madeleine (2008)
22!     dust removal, simplification by Robin Wordsworth (2009)
23!
24!     Input
25!     -----
26!     ngrid             Number of horizontal gridpoints
27!     nlayer            Number of layers
28!     nq                Number of tracers
29!     pplev             Pressure (Pa) at each layer boundary
30!     pq                Aerosol mixing ratio
31!     reffrad(ngrid,nlayer,naerkind)         Aerosol effective radius
32!     QREFvis3d(ngrid,nlayermx,naerkind) \ 3d extinction coefficients
33!     QREFir3d(ngrid,nlayermx,naerkind)  / at reference wavelengths
34!
35!     Output
36!     ------
37!     aerosol            Aerosol optical depth in layer l, grid point ig
38!     tau_col            Total column optical depth at grid point ig
39!
40!=======================================================================
41
42#include "dimensions.h"
43#include "dimphys.h"
44#include "callkeys.h"
45#include "comcstfi.h"
46#include "comvert.h"
47
48      INTEGER ngrid,nlayer,nq
49
50      REAL pplay(ngrid,nlayer)
51      REAL pplev(ngrid,nlayer+1)
52      REAL pq(ngrid,nlayer,nq)
53      REAL aerosol(ngrid,nlayer,naerkind)
54      REAL reffrad(ngrid,nlayer,naerkind)
55      REAL QREFvis3d(ngrid,nlayermx,naerkind)
56      REAL QREFir3d(ngrid,nlayermx,naerkind)
57
58      REAL tau_col(ngrid)
59!      REAL tauref(ngrid), tau_col(ngrid)
60
61      real cloudfrac(ngrid,nlayermx)
62      real aerosol0
63
64      INTEGER l,ig,iq,iaer
65
66      LOGICAL firstcall
67      DATA firstcall/.true./
68      SAVE firstcall
69
70      REAL CBRT
71      EXTERNAL CBRT
72
73      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
74      INTEGER,SAVE :: i_h2oice=0      ! water ice
75      CHARACTER(LEN=20) :: tracername ! to temporarily store text
76
77      ! for fixed dust profiles
78      real topdust, expfactor, zp
79      REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
80      REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
81      ! BENJAMIN MODIFS
82      real CLFtot
83      real totcloudfrac(ngrid)
84      logical clearsky
85
86      ! identify tracers
87      IF (firstcall) THEN
88
89        write(*,*) "Tracers found in aeropacity:"
90        do iq=1,nq
91          tracername=noms(iq)
92          if (tracername.eq."co2_ice") then
93            i_co2ice=iq
94          write(*,*) "i_co2ice=",i_co2ice
95
96          endif
97          if (tracername.eq."h2o_ice") then
98            i_h2oice=iq
99            write(*,*) "i_h2oice=",i_h2oice
100          endif
101        enddo
102
103        if (noaero) then
104          print*, "No active aerosols found in aeropacity"
105        else
106          print*, "If you would like to use aerosols, make sure any old"
107          print*, "start files are updated in newstart using the option"
108          print*, "q=0"
109          write(*,*) "Active aerosols found in aeropacity:"
110        endif
111
112        if ((iaero_co2.ne.0).and.(.not.noaero)) then
113          print*, 'iaero_co2=  ',iaero_co2
114        endif
115        if (iaero_h2o.ne.0) then
116          print*,'iaero_h2o=  ',iaero_h2o   
117        endif
118        if (iaero_dust.ne.0) then
119          print*,'iaero_dust= ',iaero_dust
120        endif
121        if (iaero_h2so4.ne.0) then
122          print*,'iaero_h2so4= ',iaero_h2so4
123        endif
124
125        firstcall=.false.
126      ENDIF ! of IF (firstcall)
127
128
129!     ---------------------------------------------------------
130!==================================================================
131!    CO2 ice aerosols
132!==================================================================
133
134      if (iaero_co2.ne.0) then
135           iaer=iaero_co2
136!       1. Initialization
137            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
138!       2. Opacity calculation
139            if (noaero) then ! aerosol set to a very low value
140             aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9
141            elseif (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
142               aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9
143               !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
144            else
145               DO ig=1, ngrid
146                  DO l=1,nlayer-1 ! to stop the rad tran bug
147
148                     aerosol0 =                         &
149                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
150                          ( rho_co2 * reffrad(ig,l,iaer) )  ) *   &
151                          ( pq(ig,l,i_co2ice) + 1.E-9 ) *         &
152                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
153                     aerosol0           = max(aerosol0,1.e-9)
154                     aerosol0           = min(aerosol0,L_TAUMAX)
155                     aerosol(ig,l,iaer) = aerosol0
156!                     aerosol(ig,l,iaer) = 0.0
157!                     print*, aerosol(ig,l,iaer)
158!        using cloud fraction
159!                     aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
160!                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
161
162
163                  ENDDO
164               ENDDO
165            end if ! if fixed or varying
166      end if ! if CO2 aerosols   
167!==================================================================
168!     Water ice / liquid
169!==================================================================
170
171      if (iaero_h2o.ne.0) then
172           iaer=iaero_h2o
173!       1. Initialization
174            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
175!       2. Opacity calculation
176            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
177               aerosol(1:ngrid,1:nlayermx,iaer) =1.e-9
178
179               ! put cloud at cloudlvl
180               if(kastprof.and.(cloudlvl.ne.0.0))then
181                  ig=1
182                  do l=1,nlayer
183                     if(int(cloudlvl).eq.l)then
184                     !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
185                        print*,'Inserting cloud at level ',l
186                        !aerosol(ig,l,iaer)=10.0
187
188                        rho_ice=920.0
189
190                        ! the Kasting approximation
191                        aerosol(ig,l,iaer) =                      &
192                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
193                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
194                          !( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
195                          ( 4.0e-4 + 1.E-9 ) *         &
196                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
197
198
199                        open(115,file='clouds.out',form='formatted')
200                        write(115,*) l,aerosol(ig,l,iaer)
201                        close(115)
202
203                        return
204                     endif
205                  end do
206
207                  call abort
208               endif
209
210            else
211
212               do ig=1, ngrid
213                  do l=1,nlayer-1 ! to stop the rad tran bug
214
215                     aerosol(ig,l,iaer) =                                    & !modification by BC
216                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
217                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
218                          !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
219                          !( pplev(ig,l) - pplev(ig,l+1) ) / g       !   opacity in the clearsky=true and the
220                                                                     !   clear=false/pq=0 case
221                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
222                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
223
224                  enddo
225               enddo
226
227               if(CLFvarying)then
228                  call totalcloudfrac(ngrid,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngrid,1:nlayermx,iaer))
229                  do ig=1, ngrid
230                     do l=1,nlayer-1 ! to stop the rad tran bug
231                        CLFtot  = max(totcloudfrac(ig),0.01)
232                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
233                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
234                     enddo
235                  enddo
236               else
237                  do ig=1, ngrid
238                     do l=1,nlayer-1 ! to stop the rad tran bug
239                        CLFtot  = CLFfixval
240                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
241                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
242                     enddo
243                  enddo
244              end if!(CLFvarying)
245            endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
246             
247      end if ! End if h2o aerosol
248
249!==================================================================
250!             Dust
251!==================================================================
252      if (iaero_dust.ne.0) then
253          iaer=iaero_dust
254!         1. Initialization
255          aerosol(1:ngrid,1:nlayermx,iaer)=0.0
256         
257          topdust=30.0 ! km  (used to be 10.0 km) LK
258
259!       2. Opacity calculation
260
261!           expfactor=0.
262           DO l=1,nlayer-1
263             DO ig=1,ngrid
264!             Typical mixing ratio profile
265
266                 zp=(preff/pplay(ig,l))**(70./topdust)
267                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
268
269!             Vertical scaling function
270              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
271               *expfactor
272
273
274             ENDDO
275           ENDDO
276
277!          Rescaling each layer to reproduce the choosen (or assimilated)
278!          dust extinction opacity at visible reference wavelength, which
279!          is scaled to the "preff" reference surface pressure available in comvert.h
280!          and stored in startfi.nc
281
282            taudusttmp(1:ngrid)=0.
283              DO l=1,nlayer
284                DO ig=1,ngrid
285                   taudusttmp(ig) = taudusttmp(ig) &
286                          +  aerosol(ig,l,iaer)
287                ENDDO
288              ENDDO
289            DO l=1,nlayer-1
290               DO ig=1,ngrid
291                  aerosol(ig,l,iaer) = max(1E-20, &
292                          dusttau &
293                       *  pplev(ig,1) / preff &
294                       *  aerosol(ig,l,iaer) &
295                       /  taudusttmp(ig))
296
297              ENDDO
298            ENDDO
299      end if ! If dust aerosol   
300
301!==================================================================
302!           H2SO4
303!==================================================================
304! added by LK
305      if (iaero_h2so4.ne.0) then
306         iaer=iaero_h2so4
307
308!       1. Initialization
309         aerosol(1:ngrid,1:nlayermx,iaer)=0.0
310
311
312!       2. Opacity calculation
313
314!           expfactor=0.
315         DO l=1,nlayer-1
316            DO ig=1,ngrid
317!              Typical mixing ratio profile
318
319               zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
320               expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
321
322!             Vertical scaling function
323               aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
324
325            ENDDO
326         ENDDO
327         tauh2so4tmp(1:ngrid)=0.
328         DO l=1,nlayer
329            DO ig=1,ngrid
330               tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
331            ENDDO
332         ENDDO
333         DO l=1,nlayer-1
334            DO ig=1,ngrid
335               aerosol(ig,l,iaer) = max(1E-20, &
336                          1 &
337                       *  pplev(ig,1) / preff &
338                       *  aerosol(ig,l,iaer) &
339                       /  tauh2so4tmp(ig))
340
341            ENDDO
342         ENDDO
343
344! 1/700. is assuming a "sulfurtau" of 1
345! Sulfur aerosol routine to be improved.
346!                     aerosol0 =                         &
347!                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
348!                          ( rho_h2so4 * reffrad(ig,l,iaer) )  ) *   &
349!                          ( pq(ig,l,i_h2so4) + 1.E-9 ) *         &
350!                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
351!                     aerosol0           = max(aerosol0,1.e-9)
352!                     aerosol0           = min(aerosol0,L_TAUMAX)
353!                     aerosol(ig,l,iaer) = aerosol0
354
355!                  ENDDO
356!               ENDDO
357      end if
358 
359           
360
361! --------------------------------------------------------------------------
362! Column integrated visible optical depth in each point (used for diagnostic)
363
364      tau_col(:)=0.0
365      do iaer = 1, naerkind
366         do l=1,nlayer
367            do ig=1,ngrid
368               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
369            end do
370         end do
371      end do
372
373      do ig=1,ngrid
374         do l=1,nlayer
375            do iaer = 1, naerkind
376               if(aerosol(ig,l,iaer).gt.1.e3)then
377                  print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
378                  print*,'at ig=',ig,',  l=',l,', iaer=',iaer
379                  print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
380                  print*,'reffrad=',reffrad(ig,l,iaer)
381               endif
382            end do
383         end do
384      end do
385
386      do ig=1,ngrid
387         if(tau_col(ig).gt.1.e3)then
388            print*,'WARNING: tau_col=',tau_col(ig)
389            print*,'at ig=',ig
390            print*,'aerosol=',aerosol(ig,:,:)
391            print*,'QREFvis3d=',QREFvis3d(ig,:,:)
392            print*,'reffrad=',reffrad(ig,:,:)
393         endif
394      end do
395      return
396    end subroutine aeropacity
397     
Note: See TracBrowser for help on using the repository browser.