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

Last change on this file since 437 was 374, checked in by emillour, 14 years ago

Generic GCM: Upgrade: The location of the 'datagcm' directory can now be given in the callphys.def file ( datadir = /absolute/path/to/datagcm ). Changed "datafile.h" into a F90 module "datafile_mod.F90" and spread this change to all routines that used to use "datafile.h".
EM

  • Property svn:executable set to *
File size: 29.7 KB
Line 
1      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,    &
2          albedo,emis,mu0,pplev,pplay,pt,                      &
3          tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
4          dtlw,dtsw,fluxsurf_lw,                               &
5          fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
6          reffrad,tau_col,cloudfrac,totcloudfrac,   &
7          clearsky,firstcall,lastcall)
8
9
10      use radinc_h
11      use radcommon_h
12      use watercommon_h
13      use datafile_mod, only: datadir
14      use ioipsl_getincom
15
16      implicit none
17
18!==================================================================
19!
20!     Purpose
21!     -------
22!     Solve the radiative transfer using the correlated-k method for
23!     the gaseous absorption and the Toon et al. (1989) method for
24!     scatttering due to aerosols.
25!
26!     Authors
27!     -------
28!     Emmanuel 01/2001, Forget 09/2001
29!     Robin Wordsworth (2009)
30!
31!==================================================================
32
33#include "dimphys.h"
34#include "comcstfi.h"
35#include "callkeys.h"
36#include "tracer.h"
37#include "gases.h"
38
39!-----------------------------------------------------------------------
40!     Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
41!     Layer #1 is the layer near the ground.
42!     Layer #nlayermx is the layer at the top.
43
44!     INPUT
45      INTEGER icount
46      INTEGER ngrid,nlayer
47      REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg)
48      REAL albedo(ngrid)                    ! SW albedo
49      REAL emis(ngrid)                      ! LW emissivity
50      REAL pplay(ngrid,nlayermx)            ! pres. level in GCM mid of layer
51      REAL pplev(ngrid,nlayermx+1)          ! pres. level at GCM layer boundaries
52
53      REAL pt(ngrid,nlayermx)               ! air temperature (K)
54      REAL tsurf(ngrid)                     ! surface temperature (K)
55      REAL dist_star,mu0(ngrid)             ! distance star-planet (AU)
56      REAL fract(ngrid)                     ! fraction of day
57
58!     Globally varying aerosol optical properties on GCM grid
59!     Not needed everywhere so not in radcommon_h
60      REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
61      REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
62      REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
63
64      REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
65      REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
66      REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
67
68      REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
69      REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
70
71!      REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
72!      REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) ! not sure of the point of these...
73
74      REAL reffrad(ngrid,nlayer,naerkind)
75      REAL nueffrad(ngrid,nlayer,naerkind)
76
77!     OUTPUT
78      REAL dtsw(ngridmx,nlayermx) ! heating rate (K/s) due to SW
79      REAL dtlw(ngridmx,nlayermx) ! heating rate (K/s) due to LW
80      REAL fluxsurf_lw(ngridmx)   ! incident LW flux to surf (W/m2)
81      REAL fluxtop_lw(ngridmx)    ! outgoing LW flux to space (W/m2)
82      REAL fluxsurf_sw(ngridmx)   ! incident SW flux to surf (W/m2)
83      REAL fluxabs_sw(ngridmx)    ! SW flux absorbed by planet (W/m2)
84      REAL fluxtop_dn(ngridmx)    ! incident top of atmosphere SW flux (W/m2)
85
86!-----------------------------------------------------------------------
87!     Declaration of the variables required by correlated-k subroutines
88!     Numbered from top to bottom unlike in the GCM!
89
90      REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
91      REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
92
93!     Optical values for the optci/cv subroutines
94      REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
95      REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
96      REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
97      REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
98      REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
99      REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
100      REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
101      REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
102      REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
103      REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
104
105      REAL*8 tauaero(L_LEVELS+1,naerkind)
106      REAL*8 nfluxtopv,nfluxtopi,nfluxtop
107      real*8 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2)
108      real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2)
109      real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
110      REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
111      REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
112      REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
113      REAL*8 albi,albv,acosz
114
115      INTEGER ig,l,k,nw,iaer,irad
116
117      real fluxtoplanet
118      real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
119      real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
120
121      real*8 qvar(L_LEVELS)          ! mixing ratio of variable component (mol/mol)
122      REAL pq(ngridmx,nlayer,nq)
123      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
124      integer nq
125
126!     Local aerosol optical properties for each column on RADIATIVE grid
127      real*8  QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
128      real*8  QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
129      real*8  GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
130      real*8  QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
131      real*8  QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
132      real*8  GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
133
134      save qxvaer, qsvaer, gvaer
135      save qxiaer, qsiaer, giaer
136      save QREFvis3d, QREFir3d
137
138      REAL tau_col(ngrid) ! diagnostic from aeropacity
139
140!     Misc.
141      logical firstcall, lastcall, nantest
142      real*8  tempv(L_NSPECTV)
143      real*8  tempi(L_NSPECTI)
144      real*8  temp,temp1,temp2,pweight
145      character(len=10) :: tmp1
146      character(len=10) :: tmp2
147
148!     for fixed water vapour profiles
149      integer i_var
150      real RH
151      real*8 pq_temp(nlayer)
152      real ptemp, Ttemp, qsat
153
154!      real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
155
156!     for OLR spec
157      integer OLRcount
158      save OLRcount
159      integer OLRcount2
160      save OLRcount2
161      character(len=2)  :: tempOLR
162      character(len=30) :: filenomOLR
163      !real ptime, pday
164      logical OLRz
165      real OLR_nu(ngrid,L_NSPECTI)
166      !real GSR_nu(ngrid,L_NSPECTV)
167      real OSR_nu(ngrid,L_NSPECTV)
168      real*8 NFLUXGNDV_nu(L_NSPECTV)
169
170!     for H2O cloud fraction in aeropacity
171      real cloudfrac(ngridmx,nlayermx)
172      real totcloudfrac(ngridmx)
173      logical clearsky
174
175!     Allow variations in cp with location
176      real cpp3D(ngridmx,nlayermx)   ! specific heat capacity at const. pressure
177
178      ! for weird cloud test
179      real pqtest(ngridmx,nlayer,nq)
180
181!     are particle radii fixed?
182      logical radfixed
183      real maxrad, minrad
184
185      real CBRT
186      external CBRT
187
188!     included by RW for runaway greenhouse 1D study
189      real muvar(ngridmx,nlayermx+1)
190      real vtmp(nlayermx)
191      REAL*8 muvarrad(L_LEVELS)
192
193      radfixed=.false.
194
195!=======================================================================
196!     Initialization on first call
197
198      qxvaer(:,:,:)=0.0
199      qsvaer(:,:,:)=0.0
200      gvaer(:,:,:) =0.0
201
202      qxiaer(:,:,:)=0.0
203      qsiaer(:,:,:)=0.0
204      giaer(:,:,:) =0.0
205
206      if(firstcall) then
207
208         call system('rm -f surf_vals_long.out')
209
210!--------------------------------------------------
211!     Effective radius and variance of the aerosols
212
213!     CO2 ice:
214         DO l=1,nlayer
215            DO ig=1,ngrid
216               reffrad(ig,l,1)  = 1.e-4
217               nueffrad(ig,l,1) = 0.1
218!     these values will change once the microphysics gets to work
219!     UNLESS tracer=.false., in which case we should be working with
220!     a fixed aerosol layer, and be able to define reffrad in a
221!     .def file. To be improved!
222            ENDDO
223         ENDDO
224
225!     H2O ice:
226         if(naerkind.eq.2)then
227            DO l=1,nlayer
228               DO ig=1,ngrid
229                  reffrad(ig,l,naerkind)  = 1.e-5
230                  nueffrad(ig,l,naerkind) = 0.1
231               ENDDO 
232            ENDDO
233         endif
234
235         print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
236         call getin("corrkdir",corrkdir)
237         print*, "corrkdir = ",corrkdir
238         write( tmp1, '(i3)' ) L_NSPECTI
239         write( tmp2, '(i3)' ) L_NSPECTV
240         banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
241         banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
242
243         call sugas_corrk       ! set up gaseous absorption properties
244         call setspi            ! basic infrared properties
245         call setspv            ! basic visible properties
246         call suaer_corrk       ! set up aerosol optical properties
247
248         Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
249
250
251         if((igcm_h2o_vap.eq.0) .and. varactive)then
252            print*,'varactive in callcorrk but no h2o_vap tracer.'
253            stop
254         endif
255
256         firstcall=.false.   
257
258      end if
259
260!=======================================================================
261!     Initialization on every call   
262
263      do l=1,nlayer
264         do ig=1,ngrid
265            do iaer=1,naerkind
266               nueffrad(ig,l,iaer) = 0.1 ! stays at 0.1
267            enddo
268         enddo
269      enddo
270
271
272      if(kastprof)then
273         radfixed=.true.
274      endif
275
276      if(radfixed)then
277         do l=1,nlayer
278            do ig=1,ngrid
279               reffrad(ig,l,1)    = 5.e-5 ! CO2 ice
280            enddo
281         enddo
282         print*,'CO2 ice particle size = ',reffrad(1,1,1)/1.e-6,' um'
283         if(naerkind.ge.2)then
284            do l=1,nlayer
285               do ig=1,ngrid
286                  !reffrad(ig,l,2) = 2.e-5 ! H2O ice
287                  reffrad(ig,l,2) = 5.e-6 ! H2O ice
288               enddo
289            enddo
290            print*,'H2O ice particle size = ',reffrad(1,1,2)/1.e-6,' um'
291         endif
292         if(naerkind.eq.3)then
293            do l=1,nlayer
294               do ig=1,ngrid
295                  reffrad(ig,l,naerkind) = 2.e-6 ! dust
296               enddo
297            enddo
298            print*,'Dust particle size = ',reffrad(1,1,naerkind)/1.e-6,' um'
299         endif
300         if(naerkind.gt.3)then
301            print*,'Code not general enough to deal with naerkind > 3 yet.'
302            call abort
303         endif
304
305      else
306
307         maxrad=0.0
308         !averad=0.0
309         minrad=1.0
310         do l=1,nlayer
311
312            !masse = (pplev(ig,l) - pplev(ig,l+1))/g
313
314            do ig=1,ngrid
315               if(tracer)then
316                  reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
317                       (4*Nmix_co2*pi*rho_co2) )
318               endif
319               reffrad(ig,l,1) = max(reffrad(ig,l,1),1.e-6)
320               reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6)
321
322               !averad = averad + reffrad(ig,l,1)*area(ig)
323               maxrad = max(reffrad(ig,l,1),maxrad)
324               minrad = min(reffrad(ig,l,1),minrad)
325            enddo
326         enddo
327         if(igcm_co2_ice.gt.0)then
328            print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
329            print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
330         endif
331
332         if((naerkind.ge.2).and.water)then
333            maxrad=0.0
334            minrad=1.0
335            do l=1,nlayer
336               do ig=1,ngrid
337                  reffrad(ig,l,2) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
338                       (4*Nmix_h2o*pi*rho_ice) )
339                  reffrad(ig,l,2) = max(reffrad(ig,l,2),1.e-6)
340                  reffrad(ig,l,2) = min(reffrad(ig,l,2),100.e-6)
341
342                  maxrad = max(reffrad(ig,l,2),maxrad)
343                  minrad = min(reffrad(ig,l,2),minrad)
344               enddo
345            enddo
346            print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
347            print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
348         endif
349
350         if(naerkind.eq.3)then
351            do l=1,nlayer
352               do ig=1,ngrid
353                  reffrad(ig,l,naerkind) = 2.e-6 ! dust
354               enddo
355            enddo
356         endif
357
358      endif
359
360
361!     how much light we get
362      fluxtoplanet=0
363      do nw=1,L_NSPECTV
364         stel(nw)=stellarf(nw)/(dist_star**2)
365         fluxtoplanet=fluxtoplanet + stel(nw)
366      end do
367
368      call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         &
369           QVISsQREF3d,omegaVIS3d,gVIS3d,                          &
370           QIRsQREF3d,omegaIR3d,gIR3d,                             &
371           QREFvis3d,QREFir3d)                                     ! get 3D aerosol optical properties
372
373      call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,      &
374           reffrad,QREFvis3d,QREFir3d,                             &
375           tau_col,cloudfrac,totcloudfrac,clearsky)                ! get aerosol optical depths
376
377
378!-----------------------------------------------------------------------
379!     Starting Big Loop over every GCM column
380      do ig=1,ngridmx
381
382!=======================================================================
383!     Transformation of the GCM variables
384
385!-----------------------------------------------------------------------
386!     Aerosol optical properties Qext, Qscat and g
387!     The transformation in the vertical is the same as for temperature
388           
389!     shortwave
390            do iaer=1,naerkind
391               DO nw=1,L_NSPECTV
392                  do l=1,nlayermx
393
394                     temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer)         &
395                         *QREFvis3d(ig,nlayermx+1-l,iaer)
396
397                     temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer)    &
398                         *QREFvis3d(ig,max(nlayermx-l,1),iaer)
399
400                     qxvaer(2*l,nw,iaer)  = temp1
401                     qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
402
403                     temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
404                     temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)
405
406                     qsvaer(2*l,nw,iaer)  = temp1
407                     qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
408
409                     temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
410                     temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)
411
412                     gvaer(2*l,nw,iaer)  = temp1
413                     gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
414
415                  end do
416
417                  qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
418                  qxvaer(2*nlayermx+1,nw,iaer)=0.
419
420                  qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
421                  qsvaer(2*nlayermx+1,nw,iaer)=0.
422
423                  gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
424                  gvaer(2*nlayermx+1,nw,iaer)=0.
425
426               end do
427
428
429!     longwave
430               DO nw=1,L_NSPECTI
431                  do l=1,nlayermx
432
433                     temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer)         &
434                          *QREFir3d(ig,nlayermx+1-l,iaer)
435
436                     temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer)    &
437                          *QREFir3d(ig,max(nlayermx-l,1),iaer)
438
439                     qxiaer(2*l,nw,iaer)  = temp1
440                     qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
441
442                     temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
443                     temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)
444
445                     qsiaer(2*l,nw,iaer)  = temp1
446                     qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
447
448                     temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
449                     temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)
450
451                     giaer(2*l,nw,iaer)  = temp1
452                     giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
453
454                  end do
455
456                  qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
457                  qxiaer(2*nlayermx+1,nw,iaer)=0.
458
459                  qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
460                  qsiaer(2*nlayermx+1,nw,iaer)=0.
461
462                  giaer(1,nw,iaer)=giaer(2,nw,iaer)
463                  giaer(2*nlayermx+1,nw,iaer)=0.
464
465               end do
466            end do
467
468
469            ! test / correct for freaky s. s. albedo values
470            do iaer=1,naerkind
471               do k=1,L_LEVELS+1
472
473                  do nw=1,L_NSPECTV
474                     if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
475                        print*,'Serious problems with qsvaer values in callcorrk'
476                        call abort
477                     endif
478                     if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
479                        qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
480                     endif
481                  end do
482
483                  do nw=1,L_NSPECTI
484                     if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
485                        print*,'Serious problems with qsiaer values in callcorrk'
486                        call abort
487                     endif
488                     if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
489                        qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
490                     endif
491                  end do
492
493               end do
494            end do
495
496!-----------------------------------------------------------------------
497!     Aerosol optical depths
498           
499         do iaer=1,naerkind     ! a bug was here           
500            do k=0,nlayer-1
501               
502               pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
503                        (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
504
505               temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
506
507               tauaero(2*k+2,iaer)=max(temp*pweight,0.0)
508               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.0)
509!
510            end do
511            ! boundary conditions
512            tauaero(1,iaer)          = tauaero(2,iaer)
513            tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
514            !tauaero(1,iaer)          = 0.
515            !tauaero(L_LEVELS+1,iaer) = 0.
516         end do
517         !print*,'Note changed tauaero BCs in callcorrk!'
518
519!     Albedo and emissivity
520         albi=1-emis(ig)        ! longwave
521         albv=albedo(ig)        ! shortwave
522
523      if(noradsurf.and.(albv.gt.0.0))then
524         print*,'For open lower boundary in callcorrk must'
525         print*,'have surface albedo set to zero!'
526         call abort
527      endif
528
529      if(ngridmx.eq.1) then       ! fixed zenith angle 'szangle' in 1D
530         acosz = cos(pi*szangle/180.0)
531         print*,'acosz=',acosz,', szangle=',szangle
532      else
533         acosz=mu0(ig)          ! cosine of sun incident angle
534      endif
535
536
537!-----------------------------------------------------------------------
538!     Water vapour (to be generalised for other gases eventually)
539     
540      if(varactive)then
541
542         i_var=igcm_h2o_vap
543         do l=1,nlayer
544            qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
545            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2   
546            ! Average approximation as for temperature...
547         end do
548         qvar(1)=qvar(2)
549         qvar(2*nlayermx+1)=qsurf(ig,i_var)
550
551      elseif(varfixed)then
552
553         do l=1,nlayermx        ! here we will assign fixed water vapour profiles globally
554            RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
555            if(RH.lt.0.0) RH=0.0
556           
557            ptemp=pplay(ig,l)
558            Ttemp=pt(ig,l)
559            call watersat(Ttemp,ptemp,qsat)
560
561            !pq_temp(l) = qsat      ! fully saturated everywhere
562            pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
563         end do
564         
565         do l=1,nlayer
566            qvar(2*l)   = pq_temp(nlayer+1-l)
567            qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
568         end do
569         qvar(1)=qvar(2)
570
571         ! Lowest layer of atmosphere
572         RH = satval * (1 - 0.02) / 0.98
573         if(RH.lt.0.0) RH=0.0
574
575         ptemp = pplev(ig,1)
576         Ttemp = tsurf(ig)
577         call watersat(Ttemp,ptemp,qsat)
578
579         !qvar(2*nlayermx+1)=qsat      ! fully saturated everywhere
580         qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
581
582      else
583         do k=1,L_LEVELS
584            qvar(k) = 1.0D-7
585         end do
586      end if
587
588      ! IMPORTANT: Now convert from kg/kg to mol/mol
589      do k=1,L_LEVELS
590         qvar(k) = qvar(k)/epsi
591      end do
592
593
594
595!-----------------------------------------------------------------------
596!     kcm mode only
597      if(kastprof)then
598
599         DO l=1,nlayer
600            muvarrad(2*l)   = mugaz
601            muvarrad(2*l+1) = mugaz
602         END DO
603
604         do k=1,L_LEVELS
605            qvar(k) = 0.0
606         end do
607         print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
608      endif
609
610
611      if(kastprof.and.(ngasmx.gt.1))then
612
613         DO l=1,nlayer
614            muvarrad(2*l)   = muvar(ig,nlayer+2-l)
615            muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
616                                muvar(ig,max(nlayer+1-l,1)))/2
617         END DO
618     
619         muvarrad(1) = muvarrad(2)
620         muvarrad(2*nlayermx+1)=muvar(ig,1)
621
622         print*,'Recalculating qvar with VARIABLE epsi for kastprof'
623         i_var=igcm_h2o_vap
624         do l=1,nlayer
625            vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
626         end do
627
628         do l=1,nlayer
629            qvar(2*l)   = vtmp(nlayer+1-l)
630            qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
631         end do
632         qvar(1)=qvar(2)
633         qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O
634
635
636      endif
637
638
639      ! Keep values inside limits for which we have radiative transfer coefficients
640      if(L_REFVAR.gt.1)then ! there was a bug here!
641         do k=1,L_LEVELS
642            if(qvar(k).lt.wrefvar(1))then
643               qvar(k)=wrefvar(1)+1.0e-8
644            elseif(qvar(k).gt.wrefvar(L_REFVAR))then
645               qvar(k)=wrefvar(L_REFVAR)-1.0e-8
646            endif
647         end do
648      endif
649
650!-----------------------------------------------------------------------
651!     Pressure and temperature
652
653      DO l=1,nlayer
654         plevrad(2*l)   = pplay(ig,nlayer+1-l)/scalep
655         plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
656         tlevrad(2*l)   = pt(ig,nlayer+1-l)
657         tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
658      END DO
659     
660      plevrad(1) = 0
661      plevrad(2) = max(pgasmin,0.0001*plevrad(3))
662
663      tlevrad(1) = tlevrad(2)
664      tlevrad(2*nlayermx+1)=tsurf(ig)
665     
666      tmid(1) = tlevrad(2)
667      tmid(2) = tlevrad(2)
668      pmid(1) = plevrad(2)
669      pmid(2) = plevrad(2)
670     
671      DO l=1,L_NLAYRAD-1
672         tmid(2*l+1) = tlevrad(2*l+1)
673         tmid(2*l+2) = tlevrad(2*l+1)
674         pmid(2*l+1) = plevrad(2*l+1)
675         pmid(2*l+2) = plevrad(2*l+1)
676      END DO
677      pmid(L_LEVELS) = plevrad(L_LEVELS)
678      tmid(L_LEVELS) = tlevrad(L_LEVELS)
679
680
681      ! test for out-of-bounds pressure
682      if(plevrad(3).lt.pgasmin)then
683         print*,'Minimum pressure is outside the radiative'
684         print*,'transfer kmatrix bounds, exiting.'
685         call abort
686      elseif(plevrad(L_LEVELS).gt.pgasmax)then
687         print*,'Maximum pressure is outside the radiative'
688         print*,'transfer kmatrix bounds, exiting.'
689         call abort
690      endif
691
692      ! test for out-of-bounds temperature
693      do k=1,L_LEVELS
694         if(tlevrad(k).lt.tgasmin)then
695            print*,'Minimum temperature is outside the radiative'
696            print*,'transfer kmatrix bounds, exiting.'
697            print*,'WARNING, OVERRIDING FOR TEST'
698            !call abort
699         elseif(tlevrad(k).gt.tgasmax)then
700            print*,'Maximum temperature is outside the radiative'
701            print*,'transfer kmatrix bounds, exiting.'
702            print*,'WARNING, OVERRIDING FOR TEST'
703            !print*, 'T=',pt
704            !call abort
705         endif
706      enddo
707
708!=======================================================================
709!     Calling the main radiative transfer subroutines
710
711
712!-----------------------------------------------------------------------
713!     Shortwave
714
715         if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
716
717            fluxtoplanet=0.
718
719            if((ngridmx.eq.1).and.(.not.(diurnal.or.tlocked)))then
720               do nw=1,L_NSPECTV
721                  stel_fract(nw)= stel(nw) * 0.25 / acosz
722                  fluxtoplanet=fluxtoplanet + stel_fract(nw)
723                                ! globally averaged = divide by 4
724                                ! but we correct for solar zenith angle
725               end do
726            else
727               do nw=1,L_NSPECTV
728                  stel_fract(nw)= stel(nw) * fract(ig)
729                  fluxtoplanet=fluxtoplanet + stel_fract(nw)
730               end do
731            endif
732
733            call optcv(dtauv,tauv,taucumv,plevrad,                 &
734                 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
735                 tmid,pmid,taugsurf,qvar,muvarrad)
736
737
738            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
739                 acosz,stel_fract,gweight,                         &
740                 nfluxtopv,nfluxoutv_nu,nfluxgndv_nu,              &
741                 !acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu,  &
742                 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
743
744         else                          ! during the night, fluxes = 0
745            nfluxtopv       = 0.0
746            nfluxoutv_nu(:) = 0.0
747            nfluxgndv_nu(:) = 0.0
748            do l=1,L_NLAYRAD
749               fmnetv(l)=0.0
750               fluxupv(l)=0.0
751               fluxdnv(l)=0.0
752            end do
753         end if
754
755
756
757
758!-----------------------------------------------------------------------
759!     Longwave
760
761         call optci(plevrad,tlevrad,dtaui,taucumi,                  &
762              qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
763              taugsurfi,qvar,muvarrad)
764
765         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
766              wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
767              fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
768
769!-----------------------------------------------------------------------
770!     Transformation of the correlated-k code outputs
771!     (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
772
773!     Flux incident at the top of the atmosphere
774         fluxtop_dn(ig)=fluxdnv(1)
775
776         fluxtop_lw(ig)  = real(nfluxtopi)
777         fluxabs_sw(ig)  = real(-nfluxtopv)
778         fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
779         fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
780
781
782         if(fluxtop_dn(ig).lt.0.0)then
783            print*,'Achtung! fluxtop_dn has lost the plot!'
784            print*,'fluxtop_dn=',fluxtop_dn(ig)
785            print*,'acosz=',acosz
786            print*,'aerosol=',aerosol(ig,:,:)
787            print*,'temp=   ',pt(ig,:)
788            print*,'pplay=  ',pplay(ig,:)
789            call abort
790         endif
791
792!     Spectral output, for exoplanet observational comparison
793         if(specOLR)then
794            do nw=1,L_NSPECTI
795               OLR_nu(ig,nw)=nfluxtopi_nu(nw)
796            end do
797            do nw=1,L_NSPECTV
798               !GSR_nu(ig,nw)=nfluxgndv_nu(nw)
799               OSR_nu(ig,nw)=nfluxoutv_nu(nw)
800            end do
801         endif
802
803!     Finally, the heating rates
804         if(nonideal)then
805
806            DO l=2,L_NLAYRAD
807               dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
808                    *g/(cpp3D(ig,L_NLAYRAD+1-l)                &
809                    *scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
810               dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
811                    *g/(cpp3D(ig,L_NLAYRAD+1-l)                &
812                    *scalep*(plevrad(2*l+1)-plevrad(2*l-1)))           
813            END DO     
814
815!     These are values at top of atmosphere
816            dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
817                 *g/(cpp3D(ig,L_NLAYRAD)*scalep*(plevrad(3)-plevrad(1)))
818            dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
819                 *g/(cpp3D(ig,L_NLAYRAD)*scalep*(plevrad(3)-plevrad(1)))
820
821         else
822
823            DO l=2,L_NLAYRAD
824               dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
825                   *g/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
826               dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
827                   *g/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
828            END DO     
829
830!     These are values at top of atmosphere
831            dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
832                *g/(cpp*scalep*(plevrad(3)-plevrad(1)))
833            dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
834                *g/(cpp*scalep*(plevrad(3)-plevrad(1)))
835
836         endif
837
838!     ---------------------------------------------------------------
839      end do                    ! end of big loop over every GCM column (ig = 1:ngrid)
840
841
842!-----------------------------------------------------------------------
843!     Additional diagnostics
844
845!     IR spectral output, for exoplanet observational comparison
846      if(specOLR)then
847        if(ngrid.ne.1)then
848          !call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu)
849          !call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu)
850        endif
851      endif
852
853      if(lastcall.and.(ngrid.eq.1))then
854
855           print*,'Saving scalar quantities in surf_vals.out...'
856           print*,'psurf = ', pplev(1,1),' Pa'
857           open(116,file='surf_vals.out')
858           write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1),         &
859                real(-nfluxtopv),real(nfluxtopi)
860           close(116)
861
862           if(specOLR)then
863               open(117,file='OLRnu.out')
864               do nw=1,L_NSPECTI
865                  write(117,*) OLR_nu(1,nw)
866               enddo
867               close(117)
868
869               open(127,file='OSRnu.out')
870               do nw=1,L_NSPECTV
871                  write(127,*) OSR_nu(1,nw)
872               enddo
873               close(127)
874           endif
875
876!     OLR vs altitude: do it as a .txt file
877           OLRz=.false.
878           if(OLRz)then
879              print*,'saving IR vertical flux for OLRz...'
880              open(118,file='OLRz_plevs.out')
881              open(119,file='OLRz.out')
882              do l=1,L_NLAYRAD
883                 write(118,*) plevrad(2*l)
884                 do nw=1,L_NSPECTI
885                     write(119,*) fluxupi_nu(l,nw)
886                  enddo
887              enddo
888              close(118)
889              close(119)
890           endif
891
892      endif
893
894    end subroutine callcorrk
Note: See TracBrowser for help on using the repository browser.