source: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90 @ 1990

Last change on this file since 1990 was 1988, checked in by jleconte, 7 years ago

28/08/2018 == JL

We now shift the radiative model top from p=0 to the middle of the last physical layer. This is done by changing pmid and plevrad in callcorrk and some corrections need to be done in gfluxv.
This seems to get rid of the aratic temperature behavior in the last two layers of the model (especially on the night side on synchronous planets).
Additional speedup corrections have been made in gfluxi that change nothing to the result.
Finally, if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv and changes in last commit.)
This has been done for water ice in aeropacity, but same correction should probably be done for other aerosol types.

  • Property svn:executable set to *
File size: 38.8 KB
Line 
1      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
2          albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
3          tsurf,fract,dist_star,aerosol,muvar,                 &
4          dtlw,dtsw,fluxsurf_lw,                               &
5          fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
6          fluxabs_sw,fluxtop_dn,                               &
7          OLR_nu,OSR_nu,                                       &
8          tau_col,cloudfrac,totcloudfrac,                      &
9          clearsky,firstcall,lastcall)
10
11      use mod_phys_lmdz_para, only : is_master
12      use radinc_h
13      use radcommon_h
14      use watercommon_h
15      use datafile_mod, only: datadir
16      use ioipsl_getin_p_mod, only: getin_p
17      use gases_h
18      use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad
19      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay, iaero_nh3, iaero_aurora
20      USE tracer_h
21      use comcstfi_mod, only: pi, mugaz, cpp
22      use callkeys_mod, only: varactive,diurnal,tracer,water,varfixed,satval,        &
23                              kastprof,strictboundcorrk,specOLR,CLFvarying
24
25      implicit none
26
27!==================================================================
28!
29!     Purpose
30!     -------
31!     Solve the radiative transfer using the correlated-k method for
32!     the gaseous absorption and the Toon et al. (1989) method for
33!     scatttering due to aerosols.
34!
35!     Authors
36!     -------
37!     Emmanuel 01/2001, Forget 09/2001
38!     Robin Wordsworth (2009)
39!
40!==================================================================
41
42!-----------------------------------------------------------------------
43!     Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
44!     Layer #1 is the layer near the ground.
45!     Layer #nlayer is the layer at the top.
46!-----------------------------------------------------------------------
47
48
49      ! INPUT
50      INTEGER,INTENT(IN) :: ngrid                  ! Number of atmospheric columns.
51      INTEGER,INTENT(IN) :: nlayer                 ! Number of atmospheric layers.
52      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)       ! Tracers (kg/kg_of_air).
53      INTEGER,INTENT(IN) :: nq                     ! Number of tracers.
54      REAL,INTENT(IN) :: qsurf(ngrid,nq)           ! Tracers on surface (kg.m-2).
55      REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV)   ! Spectral Short Wavelengths Albedo. By MT2015
56      REAL,INTENT(IN) :: emis(ngrid)               ! Long Wave emissivity.
57      REAL,INTENT(IN) :: mu0(ngrid)                ! Cosine of sun incident angle.
58      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)     ! Inter-layer pressure (Pa).
59      REAL,INTENT(IN) :: pplay(ngrid,nlayer)       ! Mid-layer pressure (Pa).
60      REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
61      REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
62      REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
63      REAL,INTENT(IN) :: dist_star                 ! Distance star-planet (AU).
64      REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
65      REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer)   ! Fraction of clouds (%).
66      logical,intent(in) :: clearsky
67      logical,intent(in) :: firstcall              ! Signals first call to physics.
68      logical,intent(in) :: lastcall               ! Signals last call to physics.
69     
70      ! OUTPUT
71      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg).
72      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer)             ! Heating rate (K/s) due to LW radiation.
73      REAL,INTENT(OUT) :: dtsw(ngrid,nlayer)             ! Heating rate (K/s) due to SW radiation.
74      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)             ! Incident LW flux to surf (W/m2).
75      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)             ! Incident SW flux to surf (W/m2)
76      REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)          ! Absorbed SW flux by the surface (W/m2). By MT2015.
77      REAL,INTENT(OUT) :: fluxtop_lw(ngrid)              ! Outgoing LW flux to space (W/m2).
78      REAL,INTENT(OUT) :: fluxabs_sw(ngrid)              ! SW flux absorbed by the planet (W/m2).
79      REAL,INTENT(OUT) :: fluxtop_dn(ngrid)              ! Incident top of atmosphere SW flux (W/m2).
80      REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)        ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
81      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
82      REAL,INTENT(OUT) :: tau_col(ngrid)                 ! Diagnostic from aeropacity.
83      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
84      REAL,INTENT(OUT) :: totcloudfrac(ngrid)            ! Column Fraction of clouds (%).
85     
86     
87     
88     
89
90      ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.   
91      REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
92      REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
93      REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
94      REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
95      REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
96      REAL :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
97
98!      REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
99!      REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
100
101      REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:)  ! aerosol effective radius (m)
102      REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
103!$OMP THREADPRIVATE(reffrad,nueffrad)
104
105!-----------------------------------------------------------------------
106!     Declaration of the variables required by correlated-k subroutines
107!     Numbered from top to bottom (unlike in the GCM)
108!-----------------------------------------------------------------------
109
110      REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
111      REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
112
113      ! Optical values for the optci/cv subroutines
114      REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
115      REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
116      REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
117      REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
118      REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
119      REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
120      REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
121      REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
122      REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
123      REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
124
125      REAL*8 tauaero(L_LEVELS,naerkind)
126      REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
127      REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
128      REAL*8 nfluxtopi_nu(L_NSPECTI)                 ! Net band-resolved IR flux at TOA (W/m2).
129      REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI)         ! For 1D diagnostic.
130      REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
131      REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
132      REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
133      REAL*8 albi,acosz
134      REAL*8 albv(L_NSPECTV)                         ! Spectral Visible Albedo.
135
136      INTEGER ig,l,k,nw,iaer
137
138      real szangle
139      logical global1d
140      save szangle,global1d
141!$OMP THREADPRIVATE(szangle,global1d)
142      real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
143      real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
144      real*8 qvar(L_LEVELS)                    ! Mixing ratio of variable component (mol/mol).
145
146      ! Local aerosol optical properties for each column on RADIATIVE grid.
147      real*8,save,allocatable ::  QXVAER(:,:,:)
148      real*8,save,allocatable ::  QSVAER(:,:,:)
149      real*8,save,allocatable ::  GVAER(:,:,:)
150      real*8,save,allocatable ::  QXIAER(:,:,:)
151      real*8,save,allocatable ::  QSIAER(:,:,:)
152      real*8,save,allocatable ::  GIAER(:,:,:)
153
154      real, dimension(:,:,:), save, allocatable :: QREFvis3d
155      real, dimension(:,:,:), save, allocatable :: QREFir3d
156!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER,QREFvis3d,QREFir3d)
157
158
159      ! Miscellaneous :
160      real*8  temp,temp1,temp2,pweight
161      character(len=10) :: tmp1
162      character(len=10) :: tmp2
163
164      ! For fixed water vapour profiles.
165      integer i_var
166      real RH
167      real*8 pq_temp(nlayer)
168! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
169      real ptemp, Ttemp, qsat
170
171      logical OLRz
172      real*8 NFLUXGNDV_nu(L_NSPECTV)
173
174      ! Included by RW for runaway greenhouse 1D study.
175      real vtmp(nlayer)
176      REAL*8 muvarrad(L_LEVELS)
177     
178      ! Included by MT for albedo calculations.     
179      REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
180      REAL*8 surface_stellar_flux   ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
181
182
183!===============================================================
184!           I.a Initialization on first call
185!===============================================================
186
187
188      if(firstcall) then
189
190        ! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
191        if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind))
192        if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind))
193        if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind))
194        if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS,L_NSPECTI,naerkind))
195        if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS,L_NSPECTI,naerkind))
196        if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS,L_NSPECTI,naerkind))
197
198         !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
199         IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
200         IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind))
201         ! Effective radius and variance of the aerosols
202         IF(.not.ALLOCATED(reffrad)) allocate(reffrad(ngrid,nlayer,naerkind))
203         IF(.not.ALLOCATED(nueffrad)) allocate(nueffrad(ngrid,nlayer,naerkind))
204
205#ifndef MESOSCALE
206         call system('rm -f surf_vals_long.out')
207#endif
208
209         if(naerkind.gt.4)then
210            print*,'Code not general enough to deal with naerkind > 4 yet.'
211            call abort
212         endif
213         call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
214         
215         
216!--------------------------------------------------
217!             Set up correlated k
218!--------------------------------------------------
219
220
221         print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
222         call getin_p("corrkdir",corrkdir)
223         print*, "corrkdir = ",corrkdir
224         write( tmp1, '(i3)' ) L_NSPECTI
225         write( tmp2, '(i3)' ) L_NSPECTV
226         banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
227         banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
228
229         call setspi            ! Basic infrared properties.
230         call setspv            ! Basic visible properties.
231         call sugas_corrk       ! Set up gaseous absorption properties.
232         call suaer_corrk       ! Set up aerosol optical properties.
233       
234
235         if((igcm_h2o_vap.eq.0) .and. varactive)then
236            print*,'varactive in callcorrk but no h2o_vap tracer.'
237            stop
238         endif
239
240         OLR_nu(:,:) = 0.
241         OSR_nu(:,:) = 0.
242
243         if (ngrid.eq.1) then
244            PRINT*, 'Simulate global averaged conditions ?'
245            global1d = .false. ! default value
246            call getin_p("global1d",global1d)
247            write(*,*) "global1d = ",global1d
248           
249            ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
250            if (global1d.and.diurnal) then
251               print*,'if global1d is true, diurnal must be set to false'
252               stop
253            endif
254
255            if (global1d) then
256               PRINT *,'Solar Zenith angle (deg.) ?'
257               PRINT *,'(assumed for averaged solar flux S/4)'
258               szangle=60.0  ! default value
259               call getin_p("szangle",szangle)
260               write(*,*) "szangle = ",szangle
261            endif
262         endif
263
264      end if ! of if (firstcall)
265
266!=======================================================================
267!          I.b  Initialization on every call   
268!=======================================================================
269 
270      qxvaer(:,:,:)=0.0
271      qsvaer(:,:,:)=0.0
272      gvaer(:,:,:) =0.0
273
274      qxiaer(:,:,:)=0.0
275      qsiaer(:,:,:)=0.0
276      giaer(:,:,:) =0.0
277
278!--------------------------------------------------
279!     Effective radius and variance of the aerosols
280!--------------------------------------------------
281
282      do iaer=1,naerkind
283
284         if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! Treat condensed co2 particles.
285            call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
286            if (is_master) then
287               print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
288               print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
289            end if
290         end if
291         
292         if ((iaer.eq.iaero_h2o).and.water) then ! Treat condensed water particles. To be generalized for other aerosols ...
293            call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
294                             reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
295            if (is_master) then
296               print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
297               print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
298            end if
299         endif
300         
301         if(iaer.eq.iaero_dust)then
302            call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
303            if (is_master) then
304               print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
305            end if
306         endif
307         
308         if(iaer.eq.iaero_h2so4)then
309            call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
310            if (is_master) then
311               print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
312            end if
313         endif
314         
315          if(iaer.eq.iaero_back2lay)then
316            call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
317         endif
318!         if(iaer.eq.iaero_nh3)then
319!           call nh3_reffrad(ngrid,nlayer,reffrad(1,1,iaero_nh3))
320!         endif
321!         if(iaer.eq.iaero_aurora)then
322!           call aurora_reffrad(ngrid,nlayer,reffrad(1,1,iaero_aurora))
323!         endif
324       
325     end do !iaer=1,naerkind.
326
327
328      ! How much light do we get ?
329      do nw=1,L_NSPECTV
330         stel(nw)=stellarf(nw)/(dist_star**2)
331      end do
332
333      ! Get 3D aerosol optical properties.
334      call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         &
335           QVISsQREF3d,omegaVIS3d,gVIS3d,                          &
336           QIRsQREF3d,omegaIR3d,gIR3d,                             &
337           QREFvis3d,QREFir3d)                                     
338
339      ! Get aerosol optical depths.
340      call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,      &
341           reffrad,QREFvis3d,QREFir3d,                             &
342           tau_col,cloudfrac,totcloudfrac,clearsky)               
343         
344
345
346!-----------------------------------------------------------------------   
347      do ig=1,ngrid ! Starting Big Loop over every GCM column
348!-----------------------------------------------------------------------
349
350
351!=======================================================================
352!              II.  Transformation of the GCM variables
353!=======================================================================
354
355
356!-----------------------------------------------------------------------
357!    Aerosol optical properties Qext, Qscat and g.
358!    The transformation in the vertical is the same as for temperature.
359!-----------------------------------------------------------------------
360           
361           
362            do iaer=1,naerkind
363               ! Shortwave.
364               do nw=1,L_NSPECTV
365               
366                  do l=1,nlayer
367
368                     temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer)         &
369                         *QREFvis3d(ig,nlayer+1-l,iaer)
370
371                     temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
372                         *QREFvis3d(ig,max(nlayer-l,1),iaer)
373
374                     qxvaer(2*l,nw,iaer)  = temp1
375                     qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
376
377                     temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
378                     temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
379
380                     qsvaer(2*l,nw,iaer)  = temp1
381                     qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
382
383                     temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
384                     temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
385
386                     gvaer(2*l,nw,iaer)  = temp1
387                     gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
388
389                  end do ! nlayer
390
391                  qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
392                  qxvaer(2*nlayer+1,nw,iaer)=0.
393
394                  qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
395                  qsvaer(2*nlayer+1,nw,iaer)=0.
396
397                  gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
398                  gvaer(2*nlayer+1,nw,iaer)=0.
399
400               end do ! L_NSPECTV
401             
402               do nw=1,L_NSPECTI
403                  ! Longwave
404                  do l=1,nlayer
405
406                     temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer)         &
407                          *QREFir3d(ig,nlayer+1-l,iaer)
408
409                     temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
410                          *QREFir3d(ig,max(nlayer-l,1),iaer)
411
412                     qxiaer(2*l,nw,iaer)  = temp1
413                     qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
414
415                     temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
416                     temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
417
418                     qsiaer(2*l,nw,iaer)  = temp1
419                     qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
420
421                     temp1=gir3d(ig,nlayer+1-l,nw,iaer)
422                     temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
423
424                     giaer(2*l,nw,iaer)  = temp1
425                     giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
426
427                  end do ! nlayer
428
429                  qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
430                  qxiaer(2*nlayer+1,nw,iaer)=0.
431
432                  qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
433                  qsiaer(2*nlayer+1,nw,iaer)=0.
434
435                  giaer(1,nw,iaer)=giaer(2,nw,iaer)
436                  giaer(2*nlayer+1,nw,iaer)=0.
437
438               end do ! L_NSPECTI
439               
440            end do ! naerkind
441
442            ! Test / Correct for freaky s. s. albedo values.
443            do iaer=1,naerkind
444               do k=1,L_LEVELS
445
446                  do nw=1,L_NSPECTV
447                     if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
448                        print*,'Serious problems with qsvaer values'
449                        print*,'in callcorrk'
450                        call abort
451                     endif
452                     if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
453                        qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
454                     endif
455                  end do
456
457                  do nw=1,L_NSPECTI
458                     if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
459                        print*,'Serious problems with qsiaer values'
460                        print*,'in callcorrk'
461                        call abort
462                     endif
463                     if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
464                        qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
465                     endif
466                  end do
467
468               end do ! L_LEVELS
469            end do ! naerkind
470
471!-----------------------------------------------------------------------
472!     Aerosol optical depths
473!-----------------------------------------------------------------------
474           
475         do iaer=1,naerkind     ! a bug was here           
476            do k=0,nlayer-1
477               
478               pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
479                       (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
480               temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
481               tauaero(2*k+2,iaer)=max(temp*pweight,0.d0)
482               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0)
483
484            end do
485            ! boundary conditions
486            tauaero(1,iaer)          = tauaero(2,iaer)
487            !tauaero(1,iaer)          = 0.
488            !JL18 at time of testing, the two above conditions gave the same results bit for bit.
489           
490         end do ! naerkind
491
492         ! Albedo and Emissivity.
493         albi=1-emis(ig)   ! Long Wave.
494         DO nw=1,L_NSPECTV ! Short Wave loop.
495            albv(nw)=albedo(ig,nw)
496         ENDDO
497
498      if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
499         acosz = cos(pi*szangle/180.0)
500         print*,'acosz=',acosz,', szangle=',szangle
501      else
502         acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
503      endif
504
505!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
506!!! Note by JL13 : In the following, some indices were changed in the interpolations,
507!!!                so that the model results are less dependent on the number of layers !
508!!!
509!!!           ---  The older versions are commented with the comment !JL13index  ---
510!!!
511!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
512
513
514!-----------------------------------------------------------------------
515!     Water vapour (to be generalised for other gases eventually ...)
516!-----------------------------------------------------------------------
517     
518      if(varactive)then
519
520         i_var=igcm_h2o_vap
521         do l=1,nlayer
522            qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
523            qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)   
524!JL13index            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2   
525!JL13index            ! Average approximation as for temperature...
526         end do
527         qvar(1)=qvar(2)
528
529      elseif(varfixed)then
530
531         do l=1,nlayer ! Here we will assign fixed water vapour profiles globally.
532            RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
533            if(RH.lt.0.0) RH=0.0
534           
535            ptemp=pplay(ig,l)
536            Ttemp=pt(ig,l)
537            call watersat(Ttemp,ptemp,qsat)
538
539            !pq_temp(l) = qsat      ! fully saturated everywhere
540            pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
541         end do
542         
543         do l=1,nlayer
544            qvar(2*l)   = pq_temp(nlayer+1-l)
545            qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
546         end do
547         
548         qvar(1)=qvar(2)
549
550         ! Lowest layer of atmosphere
551         RH = satval * (1 - 0.02) / 0.98
552         if(RH.lt.0.0) RH=0.0
553
554!         ptemp = pplev(ig,1)
555!         Ttemp = tsurf(ig)
556!         call watersat(Ttemp,ptemp,qsat)
557
558         qvar(2*nlayer+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
559 
560      else
561         do k=1,L_LEVELS
562            qvar(k) = 1.0D-7
563         end do
564      end if ! varactive/varfixed
565
566      if(.not.kastprof)then
567         ! IMPORTANT: Now convert from kg/kg to mol/mol.
568         do k=1,L_LEVELS
569            qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
570         end do
571      end if
572
573!-----------------------------------------------------------------------
574!     kcm mode only !
575!-----------------------------------------------------------------------
576
577      if(kastprof)then
578     
579         if(.not.global1d)then ! garde-fou/safeguard added by MT (to be removed in the future)
580            write(*,*) 'You have to fix mu0, '
581            write(*,*) 'the cosinus of the solar angle'
582            stop
583         endif
584         
585         ! Initial values equivalent to mugaz.
586         DO l=1,nlayer
587            muvarrad(2*l)   = mugaz
588            muvarrad(2*l+1) = mugaz
589         END DO
590
591         if(ngasmx.gt.1)then
592
593            DO l=1,nlayer
594               muvarrad(2*l)   =  muvar(ig,nlayer+2-l)
595               muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
596                                  muvar(ig,max(nlayer+1-l,1)))/2
597            END DO
598     
599            muvarrad(1) = muvarrad(2)
600            muvarrad(2*nlayer+1) = muvar(ig,1)
601
602            print*,'Recalculating qvar with VARIABLE epsi for kastprof'
603            print*,'Assumes that the variable gas is H2O!!!'
604            print*,'Assumes that there is only one tracer'
605           
606            !i_var=igcm_h2o_vap
607            i_var=1
608           
609            if(nq.gt.1)then
610               print*,'Need 1 tracer only to run kcm1d.e'
611               stop
612            endif
613           
614            do l=1,nlayer
615               vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi))
616               !vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O !JL to be changed
617            end do
618
619            do l=1,nlayer
620               qvar(2*l)   = vtmp(nlayer+1-l)
621               qvar(2*l+1) = vtmp(nlayer+1-l)
622!               qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
623            end do
624            qvar(1)=qvar(2)
625
626            print*,'Warning: reducing qvar in callcorrk.'
627            print*,'Temperature profile no longer consistent ', &
628                   'with saturated H2O. qsat=',satval
629                   
630            do k=1,L_LEVELS
631               qvar(k) = qvar(k)*satval
632            end do
633
634         endif
635      else ! if kastprof
636         DO l=1,nlayer
637            muvarrad(2*l)   = muvar(ig,nlayer+2-l)
638            muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
639         END DO
640     
641         muvarrad(1) = muvarrad(2)
642         muvarrad(2*nlayer+1)=muvar(ig,1)         
643      endif ! if kastprof
644     
645      ! Keep values inside limits for which we have radiative transfer coefficients !!!
646      if(L_REFVAR.gt.1)then ! (there was a bug here)
647         do k=1,L_LEVELS
648            if(qvar(k).lt.wrefvar(1))then
649               qvar(k)=wrefvar(1)+1.0e-8
650            elseif(qvar(k).gt.wrefvar(L_REFVAR))then
651               qvar(k)=wrefvar(L_REFVAR)-1.0e-8
652            endif
653         end do
654      endif
655
656!-----------------------------------------------------------------------
657!     Pressure and temperature
658!-----------------------------------------------------------------------
659
660      DO l=1,nlayer
661         plevrad(2*l)   = pplay(ig,nlayer+1-l)/scalep
662         plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
663         tlevrad(2*l)   = pt(ig,nlayer+1-l)
664         tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
665      END DO
666     
667      plevrad(1) = 0.
668!      plevrad(2) = 0.   !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all.
669
670      tlevrad(1) = tlevrad(2)
671      tlevrad(2*nlayer+1)=tsurf(ig)
672     
673      pmid(1) = pplay(ig,nlayer)/scalep
674      pmid(2) =  pmid(1)
675
676      tmid(1) = tlevrad(2)
677      tmid(2) = tmid(1)
678   
679      DO l=1,L_NLAYRAD-1
680         tmid(2*l+1) = tlevrad(2*l+1)
681         tmid(2*l+2) = tlevrad(2*l+1)
682         pmid(2*l+1) = plevrad(2*l+1)
683         pmid(2*l+2) = plevrad(2*l+1)
684      END DO
685      pmid(L_LEVELS) = plevrad(L_LEVELS)
686      tmid(L_LEVELS) = tlevrad(L_LEVELS)
687
688!!Alternative interpolation:
689!         pmid(3) = pmid(1)
690!         pmid(4) = pmid(1)
691!         tmid(3) = tmid(1)
692!         tmid(4) = tmid(1)
693!      DO l=2,L_NLAYRAD-1
694!         tmid(2*l+1) = tlevrad(2*l)
695!         tmid(2*l+2) = tlevrad(2*l)
696!         pmid(2*l+1) = plevrad(2*l)
697!         pmid(2*l+2) = plevrad(2*l)
698!      END DO
699!      pmid(L_LEVELS) = plevrad(L_LEVELS-1)
700!      tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
701
702      ! Test for out-of-bounds pressure.
703      if(plevrad(3).lt.pgasmin)then
704         print*,'Minimum pressure is outside the radiative'
705         print*,'transfer kmatrix bounds, exiting.'
706         call abort
707      elseif(plevrad(L_LEVELS).gt.pgasmax)then
708         print*,'Maximum pressure is outside the radiative'
709         print*,'transfer kmatrix bounds, exiting.'
710         call abort
711      endif
712
713      ! Test for out-of-bounds temperature.
714      do k=1,L_LEVELS
715         if(tlevrad(k).lt.tgasmin)then
716            print*,'Minimum temperature is outside the radiative'
717            print*,'transfer kmatrix bounds'
718            print*,"k=",k," tlevrad(k)=",tlevrad(k)
719            print*,"tgasmin=",tgasmin
720            if (strictboundcorrk) then
721              call abort
722            else
723              print*,'***********************************************'
724              print*,'we allow model to continue with tlevrad<tgasmin'
725              print*,'  ... we assume we know what you are doing ... '
726              print*,'  ... but do not let this happen too often ... '
727              print*,'***********************************************'
728              !tlevrad(k)=tgasmin ! Used in the source function !
729            endif
730         elseif(tlevrad(k).gt.tgasmax)then
731            print*,'Maximum temperature is outside the radiative'
732            print*,'transfer kmatrix bounds, exiting.'
733            print*,"k=",k," tlevrad(k)=",tlevrad(k)
734            print*,"tgasmax=",tgasmax
735            if (strictboundcorrk) then
736              call abort
737            else
738              print*,'***********************************************'
739              print*,'we allow model to continue with tlevrad<tgasmax' 
740              print*,'  ... we assume we know what you are doing ... '
741              print*,'  ... but do not let this happen too often ... '
742              print*,'***********************************************'
743              !tlevrad(k)=tgasmax ! Used in the source function !
744            endif
745         endif
746      enddo
747      do k=1,L_NLAYRAD+1
748         if(tmid(k).lt.tgasmin)then
749            print*,'Minimum temperature is outside the radiative'
750            print*,'transfer kmatrix bounds, exiting.'
751            print*,"k=",k," tmid(k)=",tmid(k)
752            print*,"tgasmin=",tgasmin
753            if (strictboundcorrk) then
754              call abort
755            else
756              print*,'***********************************************'
757              print*,'we allow model to continue but with tmid=tgasmin'
758              print*,'  ... we assume we know what you are doing ... '
759              print*,'  ... but do not let this happen too often ... '
760              print*,'***********************************************'
761              tmid(k)=tgasmin
762            endif
763         elseif(tmid(k).gt.tgasmax)then
764            print*,'Maximum temperature is outside the radiative'
765            print*,'transfer kmatrix bounds, exiting.'
766            print*,"k=",k," tmid(k)=",tmid(k)
767            print*,"tgasmax=",tgasmax
768            if (strictboundcorrk) then
769              call abort
770            else
771              print*,'***********************************************'
772              print*,'we allow model to continue but  with tmid=tgasmin'
773              print*,'  ... we assume we know what you are doing ... '
774              print*,'  ... but do not let this happen too often ... '
775              print*,'***********************************************'
776              tmid(k)=tgasmax
777            endif
778         endif
779      enddo
780
781!=======================================================================
782!          III. Calling the main radiative transfer subroutines
783!=======================================================================
784
785
786         Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
787         glat_ig=glat(ig)
788
789!-----------------------------------------------------------------------
790!        Short Wave Part
791!-----------------------------------------------------------------------
792
793         if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
794            if((ngrid.eq.1).and.(global1d))then
795               do nw=1,L_NSPECTV
796                  stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
797               end do
798            else
799               do nw=1,L_NSPECTV
800                  stel_fract(nw)= stel(nw) * fract(ig)
801               end do
802            endif
803           
804            call optcv(dtauv,tauv,taucumv,plevrad,                 &
805                 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
806                 tmid,pmid,taugsurf,qvar,muvarrad)
807
808            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
809                 acosz,stel_fract,                                 &
810                 nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,   &
811                 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
812
813         else ! During the night, fluxes = 0.
814            nfluxtopv       = 0.0d0
815            fluxtopvdn      = 0.0d0
816            nfluxoutv_nu(:) = 0.0d0
817            nfluxgndv_nu(:) = 0.0d0
818            do l=1,L_NLAYRAD
819               fmnetv(l)=0.0d0
820               fluxupv(l)=0.0d0
821               fluxdnv(l)=0.0d0
822            end do
823         end if
824
825
826         ! Equivalent Albedo Calculation (for OUTPUT). MT2015
827         if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.       
828            surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))     
829            if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
830               DO nw=1,L_NSPECTV                 
831                  albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
832               ENDDO
833               albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
834               albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
835            else
836               albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
837            endif
838         else
839            albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
840         endif
841
842
843!-----------------------------------------------------------------------
844!        Long Wave Part
845!-----------------------------------------------------------------------
846
847         call optci(plevrad,tlevrad,dtaui,taucumi,                  &
848              qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
849              taugsurfi,qvar,muvarrad)
850
851         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
852              wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu,         &
853              fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
854
855!-----------------------------------------------------------------------
856!     Transformation of the correlated-k code outputs
857!     (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
858
859!     Flux incident at the top of the atmosphere
860         fluxtop_dn(ig)=fluxtopvdn
861
862         fluxtop_lw(ig)  = real(nfluxtopi)
863         fluxabs_sw(ig)  = real(-nfluxtopv)
864         fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
865         fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
866         
867!        Flux absorbed by the surface. By MT2015.         
868         fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
869
870         if(fluxtop_dn(ig).lt.0.0)then
871            print*,'Achtung! fluxtop_dn has lost the plot!'
872            print*,'fluxtop_dn=',fluxtop_dn(ig)
873            print*,'acosz=',acosz
874            print*,'aerosol=',aerosol(ig,:,:)
875            print*,'temp=   ',pt(ig,:)
876            print*,'pplay=  ',pplay(ig,:)
877            call abort
878         endif
879
880!     Spectral output, for exoplanet observational comparison
881         if(specOLR)then
882            do nw=1,L_NSPECTI
883               OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
884            end do
885            do nw=1,L_NSPECTV
886               !GSR_nu(ig,nw)=nfluxgndv_nu(nw)
887               OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
888            end do
889         endif
890
891!     Finally, the heating rates
892
893         DO l=2,L_NLAYRAD
894            dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
895                *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
896            dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
897                *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
898         END DO     
899
900!     These are values at top of atmosphere
901         dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
902             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
903         dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
904             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2)))
905
906
907!-----------------------------------------------------------------------   
908      end do ! End of big loop over every GCM column.
909!-----------------------------------------------------------------------
910
911
912
913!-----------------------------------------------------------------------
914!     Additional diagnostics
915!-----------------------------------------------------------------------
916
917      ! IR spectral output, for exoplanet observational comparison
918      if(lastcall.and.(ngrid.eq.1))then  ! could disable the 1D output, they are in the diagfi and diagspec... JL12
919
920         print*,'Saving scalar quantities in surf_vals.out...'
921         print*,'psurf = ', pplev(1,1),' Pa'
922         open(116,file='surf_vals.out')
923         write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1),         &
924                      real(-nfluxtopv),real(nfluxtopi)
925         close(116)
926
927
928!          USEFUL COMMENT - Do Not Remove.
929!
930!           if(specOLR)then
931!               open(117,file='OLRnu.out')
932!               do nw=1,L_NSPECTI
933!                  write(117,*) OLR_nu(1,nw)
934!               enddo
935!               close(117)
936!
937!               open(127,file='OSRnu.out')
938!               do nw=1,L_NSPECTV
939!                  write(127,*) OSR_nu(1,nw)
940!               enddo
941!               close(127)
942!           endif
943
944           ! OLR vs altitude: do it as a .txt file.
945         OLRz=.false.
946         if(OLRz)then
947            print*,'saving IR vertical flux for OLRz...'
948            open(118,file='OLRz_plevs.out')
949            open(119,file='OLRz.out')
950            do l=1,L_NLAYRAD
951               write(118,*) plevrad(2*l)
952               do nw=1,L_NSPECTI
953                  write(119,*) fluxupi_nu(l,nw)
954               enddo
955            enddo
956            close(118)
957            close(119)
958         endif
959
960      endif
961
962      ! See physiq.F for explanations about CLFvarying. This is temporary.
963      if (lastcall .and. .not.CLFvarying) then
964        IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
965        IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
966!$OMP BARRIER
967!$OMP MASTER
968        IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
969        IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
970        IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
971        IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
972!$OMP END MASTER
973!$OMP BARRIER
974        IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
975        IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
976      endif
977
978
979    end subroutine callcorrk
Note: See TracBrowser for help on using the repository browser.