source: trunk/LMDZ.TITAN/libf/phytitan/get_haze_and_cloud_opacity.F90 @ 3529

Last change on this file since 3529 was 3090, checked in by slebonnois, 14 months ago

BdeBatz? : Cleans microphysics and makes few corrections for physics
  • Property svn:executable set to *
File size: 10.1 KB
Line 
1subroutine get_haze_and_cloud_opacity(FTYPE,CTYPE,m0,m3,ich,tauext,wbar,gbar)
2
3  use datafile_mod, only : datadir
4  use radinc_h, only : L_NSPECTI, L_NSPECTV
5
6  implicit none
7
8  !============================================================================
9  !
10  !     Purpose
11  !     -------
12  !     Calculates the opacity of a GCM cell from the m0, m3
13  !     and particle type information.
14  !
15  !
16  !     INPUT:
17  !     ------
18  !     FTYPE : which fractal particle type will be used in the GCM.
19  !             This is used at the first call, and defined once for all.
20  !             FTYPE = 1 : Df = 2.00
21  !             FTYPE = 2 : Df = 2.15
22  !             FTYPE = 3 : Df = 2.30
23  !             FTYPE = 4 : Df = 2.45
24  !
25  !     CTYPE : which particle type is currently called.
26  !             CTYPE = 0 : CLOUD DROPLETS
27  !             CTYPE = 1 : FRACTAL AEROSOLS (Df = 2.00)
28  !             CTYPE = 2 : FRACTAL AEROSOLS (Df = 2.15)
29  !             CTYPE = 3 : FRACTAL AEROSOLS (Df = 2.30)
30  !             CTYPE = 4 : FRACTAL AEROSOLS (Df = 2.45)
31  !             CTYPE = 5 : SPHERICAL AEROSOLS (MODE OF SMALL AEROSOLS)
32  !
33  !     m0  = 0th order moment of particule.
34  !     m3  = 3rd order moment of particule.
35  !     ich = spectral channel [1,...,46].
36  !
37  !
38  !     OUPUT:
39  !     ------
40  !     tauext = opacity.
41  !     wbar   = average single scattering albedo.
42  !     gbar   = average asymmetry parameter.
43  !
44  !
45  !     Author
46  !     ------
47  !     B. de Batz de Trenquelleon (08/2022).
48  !
49  !============================================================================
50
51
52  !-------------------------
53  ! 0. Declarations
54  !-------------------------
55
56  ! INPUTS :
57  !---------
58  integer,intent(in) :: FTYPE
59  integer,intent(in) :: CTYPE
60
61  real*8,intent(in) :: m0
62  real*8,intent(in) :: m3
63
64  integer,intent(in) :: ich
65
66  ! OUTPUTS :
67  !----------
68  real*8,intent(out) :: tauext
69  real*8,intent(out) :: wbar
70  real*8,intent(out) :: gbar
71 
72  ! PARAMETERS :
73  !-------------
74  integer :: nwl
75  integer :: nrad
76  integer :: ii,iii,kk,kkk
77
78  parameter(nwl = L_NSPECTI + L_NSPECTV, nrad = 33)
79
80  real*8 :: xwlnv(nwl)
81 
82  ! Spherical aerosols :
83  real*8, save :: rc_as(nrad)
84  real*8, save :: ssca_as(nwl,nrad),sabs_as(nwl,nrad),sext_as(nwl,nrad)
85  real*8, save :: wbar_as(nwl,nrad),gbar_as(nwl,nrad)
86  ! Fractal aerosols :
87  real*8, save :: rc_f(nrad)
88  real*8, save :: ssca_f(nwl,nrad),sabs_f(nwl,nrad),sext_f(nwl,nrad)
89  real*8, save :: wbar_f(nwl,nrad),gbar_f(nwl,nrad)
90  ! Cloud droplets :
91  real*8, save :: rc_s(nrad)
92  real*8, save :: ssca_s(nwl,nrad),sabs_s(nwl,nrad),sext_s(nwl,nrad)
93  real*8, save :: wbar_s(nwl,nrad),gbar_s(nwl,nrad)
94 
95  ! LOCAL VARIABLES FOR INTERPOLATION :
96  !------------------------------------
97  integer :: idx
98  real*8  :: rinit,rc,step,xfact
99  character(len=100) :: aer_sph_file,aer_fra_file,cloud_file
100  logical            :: file_ok
101  logical,save       :: firstcall = .true.
102
103  ! INITIALISATION :
104  !-----------------
105  tauext = 0.0
106  wbar   = 0.0
107  gbar   = 0.0
108
109  !---------------------------------------------------
110  ! 1. First call : open and read the first time
111  !---------------------------------------------------
112 
113  ! Check spectral channel :
114  if(ich.lt.1 .or. ich.gt.nwl) then
115    print*, "ERROR : ich.lt.1 .or. ich.gt.nwl in get_haze_and_cloud_opacity"
116    call abort
117  endif
118
119  if(firstcall) then
120
121    ! Spherical aerosols :
122    aer_sph_file = trim(datadir)//'/optical_tables/table_sigma_int2_sph3.00_rm50nm.dat'
123    inquire(file=trim(aer_sph_file),exist=file_ok)
124    if(.not.file_ok) then
125      write(*,*) 'The file ',TRIM(aer_sph_file),' was not found by get_haze_and_cloud_opacity.F90 !'
126      write(*,*) 'Check in ', TRIM(datadir),' that you have the one corresponding !!'
127      call abort
128    endif
129    open(unit = 27, file = aer_sph_file)
130   
131    ! Fractal aerosols :
132    if(FTYPE.eq.1) aer_fra_file = trim(datadir)//'/optical_tables/table_sigma_int2_fr2.00_rm50nm.dat'
133    if(FTYPE.eq.2) aer_fra_file = trim(datadir)//'/optical_tables/table_sigma_int2_fr2.15_rm50nm.dat'
134    if(FTYPE.eq.3) aer_fra_file = trim(datadir)//'/optical_tables/table_sigma_int2_fr2.30_rm50nm.dat'
135    if(FTYPE.eq.4) aer_fra_file = trim(datadir)//'/optical_tables/table_sigma_int2_fr2.45_rm50nm.dat'
136    inquire(file=trim(aer_fra_file),exist=file_ok)
137    if(.not.file_ok) then
138      write(*,*) 'The file ',TRIM(aer_fra_file),' was not found by get_haze_and_cloud_opacity.F90 !'
139      write(*,*) 'Check in ', TRIM(datadir),' that you have the one corresponding !!'
140      call abort
141    endif
142    open(unit = 26, file = aer_fra_file)
143   
144    ! And clouds :
145    cloud_file = trim(datadir)//'/optical_tables/table_sigma_int2_sphere.dat'
146    inquire(file=trim(cloud_file),exist=file_ok)
147    if(.not.file_ok) then
148      write(*,*) 'The file ',TRIM(cloud_file),' was not found by get_haze_and_cloud_opacity.F90 !'
149      write(*,*) 'Check in ', TRIM(datadir),' that you have the one corresponding !!'
150      call abort
151    endif
152    open(unit = 25, file = cloud_file)
153
154    do ii = 1, nwl    ! 46 loops on wlnv
155      do kk = 1, nrad ! 33 loops on rc
156
157        read(25,*) iii, kkk, xwlnv(ii), rc_s(kk), &
158                   ssca_s(ii,kk), sabs_s(ii,kk), sext_s(ii,kk), &
159                   gbar_s(ii,kk), wbar_s(ii,kk)
160        if(ii.ne.iii .or. kk.ne.kkk) then
161          write(*,*) 'lecture unit = 25'
162          write(*,*) 'ii, iii : ', ii, iii
163          write(*,*) 'kk, kkk : ', kk, kkk
164          stop
165        endif
166
167        read(26,*) iii, kkk, xwlnv(ii), rc_f(kk), &
168                   ssca_f(ii,kk), sabs_f(ii,kk), sext_f(ii,kk), &
169                   gbar_f(ii,kk), wbar_f(ii,kk)
170        if(ii.ne.iii .or. kk.ne.kkk) then
171          write(*,*) 'lecture unit = 26'
172          write(*,*) 'ii, iii : ', ii, iii
173          write(*,*) 'kk, kkk : ', kk, kkk
174          stop
175        endif
176
177        read(27,*) iii, kkk, xwlnv(ii), rc_as(kk), &
178                   ssca_as(ii,kk), sabs_as(ii,kk), sext_as(ii,kk), &
179                   gbar_as(ii,kk), wbar_as(ii,kk)
180        if(ii.ne.iii .or. kk.ne.kkk) then
181          write(*,*) 'lecture unit = 27'
182          write(*,*) 'ii, iii : ', ii, iii
183          write(*,*) 'kk, kkk : ', kk, kkk
184          stop
185       endif
186
187      enddo ! End loop on kk
188    enddo   ! End loop on ii
189
190    if(firstcall) firstcall = .false.
191
192    close(25)
193    close(26)
194    close(27)
195  endif ! End firstcall
196
197
198  !-------------------------
199  ! 2. Interpolate
200  !-------------------------
201
202  rinit = 1.e-9                ! fin = rinit*step**(33-1) = 1.e-5
203  if(CTYPE.eq.0) rinit = 1.e-7 ! fin = rinit*step**(33-1) = 1.e-3
204
205  ! No alpha(k) !
206  ! Because it is taken into account in the optical tables ...
207  rc = (m3/m0)**(1./3.)
208 
209  step = 10.**(1./8.)
210  idx = int(log(rc/rinit) / log(step)) + 1
211
212  ! Spherical (cloud) :
213  !--------------------
214  IF(CTYPE.EQ.0) THEN
215
216    if(abs(rc_s(1)-rinit)/rc_s(1).gt.1.e-5 .or. &
217       abs(rc_s(nrad)-rinit*step**32)/rc_s(nrad).gt.1.e-5) then
218        write(*,*) 'rinit          = ',rinit
219        write(*,*) 'rc_s(1)        = ',rc_s(1)
220        write(*,*) 'rinit*step**32 = ',rinit*step**32
221        write(*,*) 'rc_s(nrad)     = ',rc_s(nrad)
222        stop
223    endif
224
225    if(rc.le.rc_s(1)) then
226      wbar   = wbar_s(ich,1)
227      gbar   = gbar_s(ich,1)
228      tauext = sext_s(ich,1)*(rc/rc_s(1))**3
229    endif
230
231    if(rc.ge.rc_s(nrad)) then
232      wbar   = wbar_s(ich,nrad)
233      gbar   = gbar_s(ich,nrad)
234      tauext = sext_s(ich,nrad)*(rc/rc_s(nrad))**3.
235    endif
236
237    if(rc.gt.rc_s(1) .and. rc.lt.rc_s(nrad)) then
238      xfact  = (log(rc)-log(rc_s(idx))) / (log(rc_s(idx+1))-log(rc_s(idx)))
239      tauext = log(sext_s(ich,idx))*(1.-xfact) + log(sext_s(ich,idx+1))*xfact
240      tauext = exp(tauext)
241      wbar   = wbar_s(ich,idx)*(1.-xfact) + wbar_s(ich,idx+1)*xfact
242      gbar   = gbar_s(ich,idx)*(1.-xfact) + gbar_s(ich,idx+1)*xfact
243    endif
244
245  ENDIF ! (CTYPE.EQ.0)
246
247  ! Spherical (haze) :
248  !-------------------
249  IF(CTYPE.EQ.5) THEN
250
251    if(abs(rc_as(1)-rinit)/rc_as(1).gt.1.e-5 .or. &
252       abs(rc_as(nrad)-rinit*step**32)/rc_as(nrad).gt.1.e-5) then
253        write(*,*) 'rinit          = ',rinit
254        write(*,*) 'rc_as(1)       = ',rc_as(1)
255        write(*,*) 'rinit*step**32 = ',rinit*step**32
256        write(*,*) 'rc_as(nrad)    = ',rc_as(nrad)
257        stop
258    endif
259
260    if(rc.le.rc_as(1)) then
261      wbar   = wbar_as(ich,1)
262      gbar   = gbar_as(ich,1)
263      tauext = sext_as(ich,1)*(rc/rc_as(1))**3
264    endif
265
266    if(rc.ge.rc_as(nrad)) then
267      wbar   = wbar_as(ich,nrad)
268      gbar   = gbar_as(ich,nrad)
269      tauext = sext_as(ich,nrad)*(rc/rc_as(nrad))**3.
270    endif
271
272    if(rc.gt.rc_as(1) .and. rc.lt.rc_as(nrad)) then
273      xfact  = (log(rc)-log(rc_as(idx))) / (log(rc_as(idx+1))-log(rc_as(idx)))
274      tauext = log(sext_as(ich,idx))*(1.-xfact) + log(sext_as(ich,idx+1))*xfact
275      tauext = exp(tauext)
276      wbar   = wbar_as(ich,idx)*(1.-xfact) + wbar_as(ich,idx+1)*xfact
277      gbar   = gbar_as(ich,idx)*(1.-xfact) + gbar_as(ich,idx+1)*xfact
278    endif
279 
280  ENDIF ! (CTYPE.EQ.5)
281
282  ! Fractal (haze) :
283  !-----------------
284  IF(CTYPE.NE.0 .AND. CTYPE.NE.5) THEN
285
286    if(abs(rc_f(1)-rinit)/rc_f(1).gt.1.e-5 .or. &
287       abs(rc_f(nrad)-rinit*step**32)/rc_f(nrad).gt.1.e-5) then
288        write(*,*) 'rinit          = ',rinit
289        write(*,*) 'rc_f(1)        = ',rc_f(1)
290        write(*,*) 'rinit*step**32 = ',rinit*step**32
291        write(*,*) 'rc_f(nrad)     = ',rc_f(nrad)
292        stop
293    endif
294
295    if(rc.le.rc_f(1)) then
296      wbar   =  wbar_f(ich,1)
297      gbar   =  gbar_f(ich,1)
298      tauext = sext_f(ich,1)*(rc/rc_f(1))**3
299    endif
300
301    if(rc.ge.rc_f(nrad)) then
302      wbar   = wbar_f(ich,nrad)
303      gbar   = gbar_f(ich,nrad)
304      tauext = sext_f(ich,nrad)*(rc/rc_f(nrad))**3.
305    endif
306   
307    if(rc.gt.rc_f(1) .and. rc.lt.rc_f(nrad)) then
308      xfact  = (log(rc)-log(rc_f(idx))) / (log(rc_f(idx+1))-log(rc_f(idx)))
309      tauext = log(sext_f(ich,idx))*(1.-xfact) + log(sext_f(ich,idx+1))*xfact
310      tauext = exp(tauext)
311      wbar   = wbar_f(ich,idx)*(1.-xfact) + wbar_f(ich,idx+1)*xfact
312      gbar   = gbar_f(ich,idx)*(1.-xfact) + gbar_f(ich,idx+1)*xfact
313    endif
314
315  ENDIF ! (CTYPE.NE.0 .AND. CTYPE.NE.5)
316
317  ! Sanity check :
318  !---------------
319  IF (CTYPE.LT.0 .AND. CTYPE.GT.5) THEN
320    write(*,*) 'CTYPE < 0 ET > 5  = ',CTYPE
321    call abort
322  ENDIF
323
324  tauext = tauext * m0
325
326  return
327
328end subroutine get_haze_and_cloud_opacity
Note: See TracBrowser for help on using the repository browser.