| 1 | subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zday, & |
|---|
| 2 | albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev, & |
|---|
| 3 | pt,tsurf,fract,dist_star, & |
|---|
| 4 | dtlw,dtsw,fluxsurf_lw, & |
|---|
| 5 | fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, & |
|---|
| 6 | fluxabs_sw,fluxtop_dn, & |
|---|
| 7 | OLR_nu,OSR_nu, & |
|---|
| 8 | int_dtaui,int_dtauv,popthi,popthv,poptti,popttv, & |
|---|
| 9 | lastcall) |
|---|
| 10 | |
|---|
| 11 | use mod_phys_lmdz_para, only : is_master |
|---|
| 12 | use radinc_h |
|---|
| 13 | use radcommon_h |
|---|
| 14 | use gases_h |
|---|
| 15 | USE tracer_h |
|---|
| 16 | use callkeys_mod, only: global1d, szangle |
|---|
| 17 | use comcstfi_mod, only: pi, mugaz, cpp |
|---|
| 18 | use callkeys_mod, only: diurnal,tracer,seashaze,corrk_recombin, & |
|---|
| 19 | strictboundcorrk,specOLR,diagdtau, & |
|---|
| 20 | tplanckmin,tplanckmax,callclouds,Fcloudy |
|---|
| 21 | use geometry_mod, only: latitude |
|---|
| 22 | |
|---|
| 23 | implicit none |
|---|
| 24 | |
|---|
| 25 | !================================================================== |
|---|
| 26 | ! |
|---|
| 27 | ! Purpose |
|---|
| 28 | ! ------- |
|---|
| 29 | ! Solve the radiative transfer using the correlated-k method for |
|---|
| 30 | ! the gaseous absorption and the Toon et al. (1989) method for |
|---|
| 31 | ! scatttering due to aerosols. |
|---|
| 32 | ! |
|---|
| 33 | ! Authors |
|---|
| 34 | ! ------- |
|---|
| 35 | ! Emmanuel 01/2001, Forget 09/2001 |
|---|
| 36 | ! Robin Wordsworth (2009) |
|---|
| 37 | ! |
|---|
| 38 | ! Modified |
|---|
| 39 | ! -------- |
|---|
| 40 | ! Jan Vatant d'Ollone (2018) |
|---|
| 41 | ! --> corrk recombining case |
|---|
| 42 | ! B. de Batz de Trenquelléon (2023) |
|---|
| 43 | ! --> Titan's haze and could optics |
|---|
| 44 | ! --> clear/dark columns method |
|---|
| 45 | ! --> optical diagnostics |
|---|
| 46 | ! |
|---|
| 47 | !================================================================== |
|---|
| 48 | |
|---|
| 49 | !----------------------------------------------------------------------- |
|---|
| 50 | ! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid |
|---|
| 51 | ! Layer #1 is the layer near the ground. |
|---|
| 52 | ! Layer #nlayer is the layer at the top. |
|---|
| 53 | !----------------------------------------------------------------------- |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | ! INPUT |
|---|
| 57 | INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns. |
|---|
| 58 | INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers. |
|---|
| 59 | REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (X/kg). |
|---|
| 60 | INTEGER,INTENT(IN) :: nq ! Number of tracers. |
|---|
| 61 | REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2). |
|---|
| 62 | REAL,INTENT(IN) :: zday ! Time elapsed since Ls=0 (sols). |
|---|
| 63 | REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015 |
|---|
| 64 | REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity. |
|---|
| 65 | REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle. |
|---|
| 66 | REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa). |
|---|
| 67 | REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa). |
|---|
| 68 | REAL,INTENT(IN) :: zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries (ref : local surf). |
|---|
| 69 | REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K). |
|---|
| 70 | REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K). |
|---|
| 71 | REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day. |
|---|
| 72 | REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU). |
|---|
| 73 | logical,intent(in) :: lastcall ! Signals last call to physics. |
|---|
| 74 | |
|---|
| 75 | ! OUTPUT |
|---|
| 76 | REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation. |
|---|
| 77 | REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation. |
|---|
| 78 | REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2). |
|---|
| 79 | REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2) |
|---|
| 80 | REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015. |
|---|
| 81 | REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2). |
|---|
| 82 | REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2). |
|---|
| 83 | REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2). |
|---|
| 84 | REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1). |
|---|
| 85 | REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1). |
|---|
| 86 | REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015 |
|---|
| 87 | REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! IR optical thickness of layers within narrowbands for diags (). |
|---|
| 88 | REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! VI optical thickness of layers within narrowbands for diags (). |
|---|
| 89 | ! Diagnostics |
|---|
| 90 | REAL,INTENT(OUT) :: popthi(ngrid,nlayer,L_NSPECTI,8) ! IR optical properties of haze within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont). |
|---|
| 91 | REAL,INTENT(OUT) :: popthv(ngrid,nlayer,L_NSPECTV,8) ! VI optical properties of haze within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont). |
|---|
| 92 | REAL,INTENT(OUT) :: poptti(ngrid,nlayer,L_NSPECTI,8) ! IR optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont). |
|---|
| 93 | REAL,INTENT(OUT) :: popttv(ngrid,nlayer,L_NSPECTV,8) ! VI optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont). |
|---|
| 94 | |
|---|
| 95 | |
|---|
| 96 | |
|---|
| 97 | !----------------------------------------------------------------------- |
|---|
| 98 | ! Declaration of the variables required by correlated-k subroutines |
|---|
| 99 | ! Numbered from top to bottom (unlike in the GCM) |
|---|
| 100 | !----------------------------------------------------------------------- |
|---|
| 101 | |
|---|
| 102 | REAL*8 tmid(L_LEVELS),pmid(L_LEVELS) |
|---|
| 103 | REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS) |
|---|
| 104 | |
|---|
| 105 | ! Optical values for the optci/cv subroutines |
|---|
| 106 | REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV) |
|---|
| 107 | ! IR |
|---|
| 108 | REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 109 | REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS) |
|---|
| 110 | REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 111 | REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 112 | ! VI |
|---|
| 113 | REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 114 | REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 115 | REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS) |
|---|
| 116 | REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 117 | REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 118 | |
|---|
| 119 | ! Temporary optical values for the optci/cv subroutines (clear column) |
|---|
| 120 | REAL*8 dtaui_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 121 | REAL*8 dtauv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 122 | REAL*8 tauv_cc(L_NLEVRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 123 | REAL*8 taucumi_cc(L_LEVELS,L_NSPECTI,L_NGAUSS) |
|---|
| 124 | REAL*8 taucumv_cc(L_LEVELS,L_NSPECTV,L_NGAUSS) |
|---|
| 125 | REAL*8 wbari_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 126 | REAL*8 wbarv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 127 | REAL*8 cosbi_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 128 | REAL*8 cosbv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 129 | ! Temporary optical values for the optci/cv subroutines (dark column) |
|---|
| 130 | REAL*8 dtaui_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 131 | REAL*8 dtauv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 132 | REAL*8 tauv_dc(L_NLEVRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 133 | REAL*8 taucumi_dc(L_LEVELS,L_NSPECTI,L_NGAUSS) |
|---|
| 134 | REAL*8 taucumv_dc(L_LEVELS,L_NSPECTV,L_NGAUSS) |
|---|
| 135 | REAL*8 wbari_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 136 | REAL*8 wbarv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 137 | REAL*8 cosbi_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS) |
|---|
| 138 | REAL*8 cosbv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS) |
|---|
| 139 | |
|---|
| 140 | ! Optical diagnostics |
|---|
| 141 | ! Haze |
|---|
| 142 | REAL*8 diag_opthi(L_LEVELS,L_NSPECTI,6) |
|---|
| 143 | REAL*8 diag_opthv(L_LEVELS,L_NSPECTV,6) |
|---|
| 144 | ! Clouds |
|---|
| 145 | REAL*8 diag_optti(L_LEVELS,L_NSPECTI,6) |
|---|
| 146 | REAL*8 diag_opttv(L_LEVELS,L_NSPECTV,6) |
|---|
| 147 | ! Temporary optical diagnostics (clear column) |
|---|
| 148 | REAL*8 diag_optti_cc(L_LEVELS,L_NSPECTI,6) |
|---|
| 149 | REAL*8 diag_opttv_cc(L_LEVELS,L_NSPECTV,6) |
|---|
| 150 | ! Temporary optical diagnostics (dark column) |
|---|
| 151 | REAL*8 diag_optti_dc(L_LEVELS,L_NSPECTI,6) |
|---|
| 152 | REAL*8 diag_opttv_dc(L_LEVELS,L_NSPECTV,6) |
|---|
| 153 | |
|---|
| 154 | |
|---|
| 155 | REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn |
|---|
| 156 | REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2). |
|---|
| 157 | REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2). |
|---|
| 158 | REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic. |
|---|
| 159 | REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD) |
|---|
| 160 | REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD) |
|---|
| 161 | REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD) |
|---|
| 162 | REAL*8 albi,acosz |
|---|
| 163 | REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo. |
|---|
| 164 | |
|---|
| 165 | INTEGER ig,l,k,nw,iq,ip,ilay,it,lev2lay,cdcolumn,ng |
|---|
| 166 | |
|---|
| 167 | LOGICAL found |
|---|
| 168 | |
|---|
| 169 | real*8 taugsurf(L_NSPECTV,L_NGAUSS-1) |
|---|
| 170 | real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1) |
|---|
| 171 | |
|---|
| 172 | ! Miscellaneous : |
|---|
| 173 | character(len=100) :: message |
|---|
| 174 | character(len=10),parameter :: subname="callcorrk" |
|---|
| 175 | |
|---|
| 176 | logical OLRz |
|---|
| 177 | real*8 NFLUXGNDV_nu(L_NSPECTV) |
|---|
| 178 | |
|---|
| 179 | ! Included by MT for albedo calculations. |
|---|
| 180 | REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation. |
|---|
| 181 | REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation. |
|---|
| 182 | |
|---|
| 183 | ! For variable haze |
|---|
| 184 | REAL*8 seashazefact(L_LEVELS) |
|---|
| 185 | |
|---|
| 186 | ! For muphys optics |
|---|
| 187 | REAL*8 pqmo(ngrid,nlayer,nmicro) ! Tracers for microphysics optics (X/m2). |
|---|
| 188 | REAL*8 i2e(ngrid,nlayer) ! int 2 ext factor ( X.kg-1 -> X.m-2 for optics ) |
|---|
| 189 | |
|---|
| 190 | ! For corr-k recombining |
|---|
| 191 | REAL*8 pqr(ngrid,L_PINT,L_REFVAR) ! Tracers for corr-k recombining (mol/mol). |
|---|
| 192 | REAL*8 fact, tmin, tmax |
|---|
| 193 | |
|---|
| 194 | LOGICAL usept(L_PINT,L_NTREF) ! mask if pfref grid point will be used |
|---|
| 195 | INTEGER inflay(L_PINT) ! nearest inferior GCM layer for pfgasref grid points |
|---|
| 196 | |
|---|
| 197 | !======================================================================= |
|---|
| 198 | ! I. Initialization on every call |
|---|
| 199 | !======================================================================= |
|---|
| 200 | |
|---|
| 201 | |
|---|
| 202 | ! How much light do we get ? |
|---|
| 203 | do nw=1,L_NSPECTV |
|---|
| 204 | stel(nw)=stellarf(nw)/(dist_star**2) |
|---|
| 205 | end do |
|---|
| 206 | |
|---|
| 207 | ! Convert (microphysical) tracers for optics: X.kg-1 --> X.m-2 |
|---|
| 208 | ! NOTE: it should be moved somewhere else: calmufi performs the same kind of |
|---|
| 209 | ! computations... waste of time... |
|---|
| 210 | i2e(:,1:nlayer) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer) |
|---|
| 211 | pqmo(:,:,:) = 0.0 |
|---|
| 212 | DO iq=1,nmicro |
|---|
| 213 | pqmo(:,:,iq) = pq(:,:,iq)*i2e(:,:) |
|---|
| 214 | ENDDO |
|---|
| 215 | |
|---|
| 216 | ! Default value for fixed species for whom vmr has been taken |
|---|
| 217 | ! into account while computing high-resolution spectra |
|---|
| 218 | if (corrk_recombin) pqr(:,:,:) = 1.0 |
|---|
| 219 | |
|---|
| 220 | !----------------------------------------------------------------------- |
|---|
| 221 | do ig=1,ngrid ! Starting Big Loop over every GCM column |
|---|
| 222 | !----------------------------------------------------------------------- |
|---|
| 223 | |
|---|
| 224 | ! Recombine reference corr-k if needed |
|---|
| 225 | if (corrk_recombin) then |
|---|
| 226 | |
|---|
| 227 | ! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we |
|---|
| 228 | ! calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values |
|---|
| 229 | ! who match the atmospheric conditions ) which is then processed as a standard pre-mix in |
|---|
| 230 | ! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%. |
|---|
| 231 | |
|---|
| 232 | ! Extract tracers for variable radiative species |
|---|
| 233 | ! Also find the nearest GCM layer under each ref pressure |
|---|
| 234 | do ip=1,L_PINT |
|---|
| 235 | |
|---|
| 236 | ilay=0 |
|---|
| 237 | found = .false. |
|---|
| 238 | do l=1,nlayer |
|---|
| 239 | if ( pplay(ig,l) .gt. 10.0**(pfgasref(ip)+2.0) ) then ! pfgasref=log(p[mbar]) |
|---|
| 240 | found=.true. |
|---|
| 241 | ilay=l |
|---|
| 242 | endif |
|---|
| 243 | enddo |
|---|
| 244 | |
|---|
| 245 | if (.not. found ) then ! set to min |
|---|
| 246 | do iq=1,L_REFVAR |
|---|
| 247 | if ( radvar_mask(iq) ) then |
|---|
| 248 | pqr(ig,ip,iq) = pq(ig,1,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol |
|---|
| 249 | endif |
|---|
| 250 | enddo |
|---|
| 251 | else |
|---|
| 252 | if (ilay==nlayer) then ! set to max |
|---|
| 253 | do iq=1,L_REFVAR |
|---|
| 254 | if ( radvar_mask(iq) ) then |
|---|
| 255 | pqr(ig,ip,iq) = pq(ig,nlayer,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol |
|---|
| 256 | endif |
|---|
| 257 | enddo |
|---|
| 258 | else ! standard |
|---|
| 259 | fact = ( 10.0**(pfgasref(ip)+2.0) - pplay(ig,ilay+1) ) / ( pplay(ig,ilay) - pplay(ig,ilay+1) ) ! pfgasref=log(p[mbar]) |
|---|
| 260 | do iq=1,L_REFVAR |
|---|
| 261 | if ( radvar_mask(iq) ) then |
|---|
| 262 | pqr(ig,ip,iq) = pq(ig,ilay,radvar_indx(iq))**fact * pq(ig,ilay+1,radvar_indx(iq))**(1.0-fact) |
|---|
| 263 | pqr(ig,ip,iq) = pqr(ig,ip,iq) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol |
|---|
| 264 | endif |
|---|
| 265 | enddo |
|---|
| 266 | endif ! if ilay==nlayer |
|---|
| 267 | endif ! if not found |
|---|
| 268 | |
|---|
| 269 | inflay(ip) = ilay |
|---|
| 270 | |
|---|
| 271 | enddo ! ip=1,L_PINT |
|---|
| 272 | |
|---|
| 273 | ! NB : The following usept is a trick to call recombine only for the reference T-P |
|---|
| 274 | ! grid points that are useful given the temperature range at this altitude |
|---|
| 275 | ! It saves a looot of time - JVO 18 |
|---|
| 276 | usept(:,:) = .true. |
|---|
| 277 | do ip=1,L_PINT-1 |
|---|
| 278 | if ( inflay(ip+1)==nlayer ) then |
|---|
| 279 | usept(ip,:) = .false. |
|---|
| 280 | endif |
|---|
| 281 | if ( inflay(ip) == 0 ) then |
|---|
| 282 | usept(ip+1:,:) = .false. |
|---|
| 283 | endif |
|---|
| 284 | if ( usept(ip,1) ) then ! if not all false |
|---|
| 285 | tmin = minval(pt(ig,max(inflay(ip+1)+1,1):min(inflay(ip)+1,nlayer))) |
|---|
| 286 | tmax = maxval(pt(ig,max(inflay(ip+1)+1,1):min(inflay(ip)+1,nlayer))) |
|---|
| 287 | do it=1,L_NTREF-1 |
|---|
| 288 | if ( tgasref(it+1) .lt. tmin ) then |
|---|
| 289 | usept(ip,it) = .false. |
|---|
| 290 | endif |
|---|
| 291 | enddo |
|---|
| 292 | do it=2,L_NTREF |
|---|
| 293 | if ( tgasref(it-1) .gt. tmax ) then |
|---|
| 294 | usept(ip,it) = .false. |
|---|
| 295 | endif |
|---|
| 296 | enddo |
|---|
| 297 | ! in case of out-of-bounds |
|---|
| 298 | if ( tgasref(1) .lt. tmin ) usept(ip,1) = .true. |
|---|
| 299 | if ( tgasref(L_NTREF) .gt. tmax ) usept(ip,L_NTREF) = .true. |
|---|
| 300 | endif |
|---|
| 301 | enddo ! ip=1,L_PINT-1 |
|---|
| 302 | ! deal with last bound |
|---|
| 303 | if ( inflay(L_PINT-1).ne.0 ) usept(L_PINT,:) = usept(L_PINT-1,:) |
|---|
| 304 | |
|---|
| 305 | |
|---|
| 306 | do ip=1,L_PINT |
|---|
| 307 | |
|---|
| 308 | ! Recombine k at (useful only!) reference T-P values if tracers or T have enough varied |
|---|
| 309 | do it=1,L_NTREF |
|---|
| 310 | |
|---|
| 311 | if ( usept(ip,it) .eqv. .false. ) cycle |
|---|
| 312 | |
|---|
| 313 | do l=1,L_REFVAR |
|---|
| 314 | if ( abs( (pqr(ig,ip,l) - pqrold(ip,l)) / max(1.0e-30,pqrold(ip,l))) .GT. 0.01 & ! +- 1% |
|---|
| 315 | .or. ( useptold(ip,it) .eqv. .false. ) ) then ! in case T change but not the tracers |
|---|
| 316 | call recombin_corrk( pqr(ig,ip,:),ip,it ) |
|---|
| 317 | exit ! one is enough to trigger the update |
|---|
| 318 | endif |
|---|
| 319 | enddo |
|---|
| 320 | |
|---|
| 321 | enddo |
|---|
| 322 | |
|---|
| 323 | enddo ! ip=1,L_PINT |
|---|
| 324 | |
|---|
| 325 | useptold(:,:)=usept(:,:) |
|---|
| 326 | |
|---|
| 327 | endif ! if corrk_recombin |
|---|
| 328 | |
|---|
| 329 | !======================================================================= |
|---|
| 330 | ! II. Transformation of the GCM variables |
|---|
| 331 | !======================================================================= |
|---|
| 332 | |
|---|
| 333 | |
|---|
| 334 | ! Albedo and Emissivity. |
|---|
| 335 | albi=1-emis(ig) ! Long Wave. |
|---|
| 336 | DO nw=1,L_NSPECTV ! Short Wave loop. |
|---|
| 337 | albv(nw)=albedo(ig,nw) |
|---|
| 338 | ENDDO |
|---|
| 339 | |
|---|
| 340 | if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight. |
|---|
| 341 | acosz = cos(pi*szangle/180.0) |
|---|
| 342 | print*,'acosz=',acosz,', szangle=',szangle |
|---|
| 343 | else |
|---|
| 344 | acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude. |
|---|
| 345 | endif |
|---|
| 346 | |
|---|
| 347 | !----------------------------------------------------------------------- |
|---|
| 348 | ! Pressure and temperature |
|---|
| 349 | !----------------------------------------------------------------------- |
|---|
| 350 | |
|---|
| 351 | DO l=1,nlayer |
|---|
| 352 | plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep |
|---|
| 353 | plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep |
|---|
| 354 | tlevrad(2*l) = pt(ig,nlayer+1-l) |
|---|
| 355 | tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2 |
|---|
| 356 | END DO |
|---|
| 357 | |
|---|
| 358 | plevrad(1) = 0. |
|---|
| 359 | plevrad(2) = 0. !! Trick to have correct calculations of fluxes in gflux(i/v).F, but the pmid levels are not impacted by this change. |
|---|
| 360 | |
|---|
| 361 | tlevrad(1) = tlevrad(2) |
|---|
| 362 | tlevrad(2*nlayer+1)=tsurf(ig) |
|---|
| 363 | |
|---|
| 364 | pmid(1) = max(pgasmin,0.0001*plevrad(3)) |
|---|
| 365 | pmid(2) = pmid(1) |
|---|
| 366 | |
|---|
| 367 | tmid(1) = tlevrad(2) |
|---|
| 368 | tmid(2) = tmid(1) |
|---|
| 369 | |
|---|
| 370 | DO l=1,L_NLAYRAD-1 |
|---|
| 371 | tmid(2*l+1) = tlevrad(2*l+1) |
|---|
| 372 | tmid(2*l+2) = tlevrad(2*l+1) |
|---|
| 373 | pmid(2*l+1) = plevrad(2*l+1) |
|---|
| 374 | pmid(2*l+2) = plevrad(2*l+1) |
|---|
| 375 | END DO |
|---|
| 376 | pmid(L_LEVELS) = plevrad(L_LEVELS) |
|---|
| 377 | tmid(L_LEVELS) = tlevrad(L_LEVELS) |
|---|
| 378 | |
|---|
| 379 | !!Alternative interpolation: |
|---|
| 380 | ! pmid(3) = pmid(1) |
|---|
| 381 | ! pmid(4) = pmid(1) |
|---|
| 382 | ! tmid(3) = tmid(1) |
|---|
| 383 | ! tmid(4) = tmid(1) |
|---|
| 384 | ! DO l=2,L_NLAYRAD-1 |
|---|
| 385 | ! tmid(2*l+1) = tlevrad(2*l) |
|---|
| 386 | ! tmid(2*l+2) = tlevrad(2*l) |
|---|
| 387 | ! pmid(2*l+1) = plevrad(2*l) |
|---|
| 388 | ! pmid(2*l+2) = plevrad(2*l) |
|---|
| 389 | ! END DO |
|---|
| 390 | ! pmid(L_LEVELS) = plevrad(L_LEVELS-1) |
|---|
| 391 | ! tmid(L_LEVELS) = tlevrad(L_LEVELS-1) |
|---|
| 392 | |
|---|
| 393 | ! Test for out-of-bounds pressure. |
|---|
| 394 | if(plevrad(3).lt.pgasmin)then |
|---|
| 395 | print*,'Minimum pressure is outside the radiative' |
|---|
| 396 | print*,'transfer kmatrix bounds, exiting.' |
|---|
| 397 | call abort |
|---|
| 398 | elseif(plevrad(L_LEVELS).gt.pgasmax)then |
|---|
| 399 | print*,'Maximum pressure is outside the radiative' |
|---|
| 400 | print*,'transfer kmatrix bounds, exiting.' |
|---|
| 401 | call abort |
|---|
| 402 | endif |
|---|
| 403 | |
|---|
| 404 | ! Test for out-of-bounds temperature. |
|---|
| 405 | do k=1,L_LEVELS |
|---|
| 406 | if(tlevrad(k).lt.tgasmin)then |
|---|
| 407 | print*,'Minimum temperature is outside the radiative' |
|---|
| 408 | print*,'transfer kmatrix bounds' |
|---|
| 409 | print*,"k=",k," tlevrad(k)=",tlevrad(k) |
|---|
| 410 | print*,"tgasmin=",tgasmin |
|---|
| 411 | if (strictboundcorrk) then |
|---|
| 412 | call abort |
|---|
| 413 | else |
|---|
| 414 | print*,'***********************************************' |
|---|
| 415 | print*,'we allow model to continue with tlevrad=tgasmin' |
|---|
| 416 | print*,' ... we assume we know what you are doing ... ' |
|---|
| 417 | print*,' ... but do not let this happen too often ... ' |
|---|
| 418 | print*,'***********************************************' |
|---|
| 419 | !tlevrad(k)=tgasmin |
|---|
| 420 | endif |
|---|
| 421 | elseif(tlevrad(k).gt.tgasmax)then |
|---|
| 422 | ! print*,'Maximum temperature is outside the radiative' |
|---|
| 423 | ! print*,'transfer kmatrix bounds, exiting.' |
|---|
| 424 | ! print*,"k=",k," tlevrad(k)=",tlevrad(k) |
|---|
| 425 | ! print*,"tgasmax=",tgasmax |
|---|
| 426 | if (strictboundcorrk) then |
|---|
| 427 | call abort |
|---|
| 428 | else |
|---|
| 429 | ! print*,'***********************************************' |
|---|
| 430 | ! print*,'we allow model to continue with tlevrad=tgasmax' |
|---|
| 431 | ! print*,' ... we assume we know what you are doing ... ' |
|---|
| 432 | ! print*,' ... but do not let this happen too often ... ' |
|---|
| 433 | ! print*,'***********************************************' |
|---|
| 434 | !tlevrad(k)=tgasmax |
|---|
| 435 | endif |
|---|
| 436 | endif |
|---|
| 437 | |
|---|
| 438 | if (tlevrad(k).lt.tplanckmin) then |
|---|
| 439 | print*,'Minimum temperature is outside the boundaries for' |
|---|
| 440 | print*,'Planck function integration set in callphys.def, aborting.' |
|---|
| 441 | print*,"k=",k," tlevrad(k)=",tlevrad(k) |
|---|
| 442 | print*,"tplanckmin=",tplanckmin |
|---|
| 443 | message="Minimum temperature outside Planck function bounds - Change tplanckmin in callphys.def" |
|---|
| 444 | call abort_physic(subname,message,1) |
|---|
| 445 | else if (tlevrad(k).gt.tplanckmax) then |
|---|
| 446 | print*,'Maximum temperature is outside the boundaries for' |
|---|
| 447 | print*,'Planck function integration set in callphys.def, aborting.' |
|---|
| 448 | print*,"k=",k," tlevrad(k)=",tlevrad(k) |
|---|
| 449 | print*,"tplanckmax=",tplanckmax |
|---|
| 450 | message="Maximum temperature outside Planck function bounds - Change tplanckmax in callphys.def" |
|---|
| 451 | call abort_physic(subname,message,1) |
|---|
| 452 | endif |
|---|
| 453 | |
|---|
| 454 | enddo |
|---|
| 455 | |
|---|
| 456 | do k=1,L_NLAYRAD+1 |
|---|
| 457 | if(tmid(k).lt.tgasmin)then |
|---|
| 458 | print*,'Minimum temperature is outside the radiative' |
|---|
| 459 | print*,'transfer kmatrix bounds, exiting.' |
|---|
| 460 | print*,"k=",k," tmid(k)=",tmid(k) |
|---|
| 461 | print*,"tgasmin=",tgasmin |
|---|
| 462 | if (strictboundcorrk) then |
|---|
| 463 | call abort |
|---|
| 464 | else |
|---|
| 465 | print*,'***********************************************' |
|---|
| 466 | print*,'we allow model to continue with tmid=tgasmin' |
|---|
| 467 | print*,' ... we assume we know what you are doing ... ' |
|---|
| 468 | print*,' ... but do not let this happen too often ... ' |
|---|
| 469 | print*,'***********************************************' |
|---|
| 470 | tmid(k)=tgasmin |
|---|
| 471 | endif |
|---|
| 472 | elseif(tmid(k).gt.tgasmax)then |
|---|
| 473 | ! print*,'Maximum temperature is outside the radiative' |
|---|
| 474 | ! print*,'transfer kmatrix bounds, exiting.' |
|---|
| 475 | ! print*,"k=",k," tmid(k)=",tmid(k) |
|---|
| 476 | ! print*,"tgasmax=",tgasmax |
|---|
| 477 | if (strictboundcorrk) then |
|---|
| 478 | call abort |
|---|
| 479 | else |
|---|
| 480 | ! print*,'***********************************************' |
|---|
| 481 | ! print*,'we allow model to continue with tmid=tgasmin' |
|---|
| 482 | ! print*,' ... we assume we know what you are doing ... ' |
|---|
| 483 | ! print*,' ... but do not let this happen too often ... ' |
|---|
| 484 | ! print*,'***********************************************' |
|---|
| 485 | tmid(k)=tgasmax |
|---|
| 486 | endif |
|---|
| 487 | endif |
|---|
| 488 | enddo |
|---|
| 489 | |
|---|
| 490 | !======================================================================= |
|---|
| 491 | ! III. Calling the main radiative transfer subroutines |
|---|
| 492 | !======================================================================= |
|---|
| 493 | |
|---|
| 494 | Cmk(:) = 0.01 * 1.0 / (gzlat(ig,:) * mugaz * 1.672621e-27) ! q_main=1.0 assumed. |
|---|
| 495 | gzlat_ig(:) = gzlat(ig,:) |
|---|
| 496 | |
|---|
| 497 | ! Compute the haze seasonal modulation factor |
|---|
| 498 | if (seashaze) call season_haze(zday,latitude(ig),plevrad,seashazefact) |
|---|
| 499 | |
|---|
| 500 | !----------------------------------------------------------------------- |
|---|
| 501 | ! Short Wave Part |
|---|
| 502 | !----------------------------------------------------------------------- |
|---|
| 503 | |
|---|
| 504 | ! Clear column : |
|---|
| 505 | cdcolumn = 0 |
|---|
| 506 | call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid, & |
|---|
| 507 | dtauv_cc,tauv_cc,taucumv_cc,wbarv_cc,cosbv_cc,tauray,taugsurf,seashazefact,& |
|---|
| 508 | diag_opthv,diag_opttv_cc,cdcolumn) |
|---|
| 509 | ! Dark column : |
|---|
| 510 | cdcolumn = 1 |
|---|
| 511 | call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid, & |
|---|
| 512 | dtauv_dc,tauv_dc,taucumv_dc,wbarv_dc,cosbv_dc,tauray,taugsurf,seashazefact,& |
|---|
| 513 | diag_opthv,diag_opttv_dc,cdcolumn) |
|---|
| 514 | |
|---|
| 515 | ! Mean opacity, ssa and asf : |
|---|
| 516 | where ((exp(-dtauv_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtauv_dc(:,:,:)).ge.1.d-40)) |
|---|
| 517 | dtauv(:,:,:) = -log((1-Fcloudy)*exp(-dtauv_cc(:,:,:)) + Fcloudy*exp(-dtauv_dc(:,:,:))) |
|---|
| 518 | elsewhere |
|---|
| 519 | dtauv(:,:,:) = dtauv_dc(:,:,:) ! No need to average... |
|---|
| 520 | endwhere |
|---|
| 521 | do ng = 1, L_NGAUSS |
|---|
| 522 | do nw = 1, L_NSPECTV |
|---|
| 523 | taucumv(1,nw,ng) = 0.0d0 |
|---|
| 524 | do k = 2, L_LEVELS |
|---|
| 525 | if ((exp(-taucumv_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumv_dc(k,nw,ng)).ge.1.d-40)) then |
|---|
| 526 | taucumv(k,nw,ng) = taucumv(k-1,nw,ng) - log((1-Fcloudy)*exp(-taucumv_cc(k,nw,ng)) + Fcloudy*exp(-taucumv_dc(k,nw,ng))) |
|---|
| 527 | else |
|---|
| 528 | taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + taucumv_dc(k,nw,ng) ! No need to average... |
|---|
| 529 | end if |
|---|
| 530 | end do |
|---|
| 531 | do l = 1, L_NLAYRAD |
|---|
| 532 | tauv(l,nw,ng) = taucumv(2*l,nw,ng) |
|---|
| 533 | end do |
|---|
| 534 | tauv(l,nw,ng) = taucumv(2*L_NLAYRAD+1,nw,ng) |
|---|
| 535 | end do |
|---|
| 536 | end do |
|---|
| 537 | wbarv = (1-Fcloudy) * wbarv_cc + Fcloudy * wbarv_dc |
|---|
| 538 | cosbv = (1-Fcloudy) * cosbv_cc + Fcloudy * cosbv_dc |
|---|
| 539 | |
|---|
| 540 | ! Diagnostics for clouds : |
|---|
| 541 | if (callclouds) then |
|---|
| 542 | where (diag_opttv_cc(:,:,1) .lt. 1.d-30) |
|---|
| 543 | diag_opttv_cc(:,:,1) = 1.d-30 |
|---|
| 544 | endwhere |
|---|
| 545 | where (diag_opttv_dc(:,:,1) .lt. 1.d-30) |
|---|
| 546 | diag_opttv_dc(:,:,1) = 1.d-30 |
|---|
| 547 | endwhere |
|---|
| 548 | diag_opttv(:,:,1) = -log( (1-Fcloudy)*exp(-diag_opttv_cc(:,:,1)) + Fcloudy*exp(-diag_opttv_dc(:,:,1)) ) |
|---|
| 549 | diag_opttv(:,:,2) = (1-Fcloudy) * diag_opttv_cc(:,:,2) + Fcloudy * diag_opttv_dc(:,:,2) |
|---|
| 550 | diag_opttv(:,:,3) = (1-Fcloudy) * diag_opttv_cc(:,:,3) + Fcloudy * diag_opttv_dc(:,:,3) |
|---|
| 551 | diag_opttv(:,:,4) = -log( (1-Fcloudy)*exp(-diag_opttv_cc(:,:,4)) + Fcloudy*exp(-diag_opttv_dc(:,:,4)) ) |
|---|
| 552 | diag_opttv(:,:,5) = -log( (1-Fcloudy)*exp(-diag_opttv_cc(:,:,5)) + Fcloudy*exp(-diag_opttv_dc(:,:,5)) ) |
|---|
| 553 | diag_opttv(:,:,6) = -log( (1-Fcloudy)*exp(-diag_opttv_cc(:,:,6)) + Fcloudy*exp(-diag_opttv_dc(:,:,6)) ) |
|---|
| 554 | endif |
|---|
| 555 | |
|---|
| 556 | if(fract(ig) .ge. 1.0e-4) then ! Only during daylight. |
|---|
| 557 | if((ngrid.eq.1).and.(global1d))then |
|---|
| 558 | do nw=1,L_NSPECTV |
|---|
| 559 | stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle |
|---|
| 560 | end do |
|---|
| 561 | else |
|---|
| 562 | do nw=1,L_NSPECTV |
|---|
| 563 | stel_fract(nw)= stel(nw) * fract(ig) |
|---|
| 564 | end do |
|---|
| 565 | endif |
|---|
| 566 | call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, & |
|---|
| 567 | acosz,stel_fract, & |
|---|
| 568 | nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, & |
|---|
| 569 | fmnetv,fluxupv,fluxdnv,fzerov,taugsurf) |
|---|
| 570 | |
|---|
| 571 | else ! During the night, fluxes = 0. |
|---|
| 572 | nfluxtopv = 0.0d0 |
|---|
| 573 | fluxtopvdn = 0.0d0 |
|---|
| 574 | nfluxoutv_nu(:) = 0.0d0 |
|---|
| 575 | nfluxgndv_nu(:) = 0.0d0 |
|---|
| 576 | do l=1,L_NLAYRAD |
|---|
| 577 | fmnetv(l)=0.0d0 |
|---|
| 578 | fluxupv(l)=0.0d0 |
|---|
| 579 | fluxdnv(l)=0.0d0 |
|---|
| 580 | end do |
|---|
| 581 | end if |
|---|
| 582 | |
|---|
| 583 | |
|---|
| 584 | ! Equivalent Albedo Calculation (for OUTPUT). MT2015 |
|---|
| 585 | if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight. |
|---|
| 586 | surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV)) |
|---|
| 587 | if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive. |
|---|
| 588 | DO nw=1,L_NSPECTV |
|---|
| 589 | albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw) |
|---|
| 590 | ENDDO |
|---|
| 591 | albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux |
|---|
| 592 | albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV)) |
|---|
| 593 | else |
|---|
| 594 | albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0. |
|---|
| 595 | endif |
|---|
| 596 | else |
|---|
| 597 | albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0. |
|---|
| 598 | endif |
|---|
| 599 | |
|---|
| 600 | |
|---|
| 601 | !----------------------------------------------------------------------- |
|---|
| 602 | ! Long Wave Part |
|---|
| 603 | !----------------------------------------------------------------------- |
|---|
| 604 | |
|---|
| 605 | |
|---|
| 606 | ! Clear column : |
|---|
| 607 | cdcolumn = 0 |
|---|
| 608 | call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,& |
|---|
| 609 | dtaui_cc,taucumi_cc,cosbi_cc,wbari_cc,taugsurfi,seashazefact, & |
|---|
| 610 | diag_opthi,diag_optti_cc,cdcolumn) |
|---|
| 611 | ! Dark column : |
|---|
| 612 | cdcolumn = 1 |
|---|
| 613 | call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,& |
|---|
| 614 | dtaui_dc,taucumi_dc,cosbi_dc,wbari_dc,taugsurfi,seashazefact, & |
|---|
| 615 | diag_opthi,diag_optti_dc,cdcolumn) |
|---|
| 616 | |
|---|
| 617 | ! Mean opacity, ssa and asf : |
|---|
| 618 | where ((exp(-dtaui_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtaui_dc(:,:,:)).ge.1.d-40)) |
|---|
| 619 | dtaui(:,:,:) = -log((1-Fcloudy)*exp(-dtaui_cc(:,:,:)) + Fcloudy*exp(-dtaui_dc(:,:,:))) |
|---|
| 620 | elsewhere |
|---|
| 621 | dtaui(:,:,:) = dtaui_dc(:,:,:) ! No need to average... |
|---|
| 622 | endwhere |
|---|
| 623 | do ng = 1, L_NGAUSS |
|---|
| 624 | do nw = 1, L_NSPECTI |
|---|
| 625 | taucumi(1,nw,ng) = 0.0d0 |
|---|
| 626 | do k = 2, L_LEVELS |
|---|
| 627 | if ((exp(-taucumi_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumi_dc(k,nw,ng)).ge.1.d-40)) then |
|---|
| 628 | taucumi(k,nw,ng) = taucumi(k-1,nw,ng) - log((1-Fcloudy)*exp(-taucumi_cc(k,nw,ng)) + Fcloudy*exp(-taucumi_dc(k,nw,ng))) |
|---|
| 629 | else |
|---|
| 630 | taucumi(k,nw,ng) = taucumi(k-1,nw,ng) + taucumi_dc(k,nw,ng) ! No need to average... |
|---|
| 631 | end if |
|---|
| 632 | end do |
|---|
| 633 | end do |
|---|
| 634 | end do |
|---|
| 635 | wbari = (1-Fcloudy) * wbari_cc + Fcloudy * wbari_dc |
|---|
| 636 | cosbi = (1-Fcloudy) * cosbi_cc + Fcloudy * cosbi_dc |
|---|
| 637 | |
|---|
| 638 | ! Diagnostics for clouds : |
|---|
| 639 | if (callclouds) then |
|---|
| 640 | where (diag_optti_cc(:,:,1) .lt. 1.d-30) |
|---|
| 641 | diag_optti_cc(:,:,1) = 1.d-30 |
|---|
| 642 | endwhere |
|---|
| 643 | where (diag_optti_dc(:,:,1) .lt. 1.d-30) |
|---|
| 644 | diag_optti_dc(:,:,1) = 1.d-30 |
|---|
| 645 | endwhere |
|---|
| 646 | diag_optti(:,:,1) = -log( (1-Fcloudy)*exp(-diag_optti_cc(:,:,1)) + Fcloudy*exp(-diag_optti_dc(:,:,1)) ) |
|---|
| 647 | diag_optti(:,:,2) = (1-Fcloudy) * diag_optti_cc(:,:,2) + Fcloudy * diag_optti_dc(:,:,2) |
|---|
| 648 | diag_optti(:,:,3) = (1-Fcloudy) * diag_optti_cc(:,:,3) + Fcloudy * diag_optti_dc(:,:,3) |
|---|
| 649 | diag_optti(:,:,4) = -log( (1-Fcloudy)*exp(-diag_optti_cc(:,:,4)) + Fcloudy*exp(-diag_optti_dc(:,:,4)) ) |
|---|
| 650 | diag_optti(:,:,5) = -log( (1-Fcloudy)*exp(-diag_optti_cc(:,:,5)) + Fcloudy*exp(-diag_optti_dc(:,:,5)) ) |
|---|
| 651 | diag_optti(:,:,6) = -log( (1-Fcloudy)*exp(-diag_optti_cc(:,:,6)) + Fcloudy*exp(-diag_optti_dc(:,:,6)) ) |
|---|
| 652 | endif |
|---|
| 653 | |
|---|
| 654 | call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, & |
|---|
| 655 | wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, & |
|---|
| 656 | fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi) |
|---|
| 657 | |
|---|
| 658 | !----------------------------------------------------------------------- |
|---|
| 659 | ! Transformation of the correlated-k code outputs |
|---|
| 660 | ! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw) |
|---|
| 661 | |
|---|
| 662 | ! Flux incident at the top of the atmosphere |
|---|
| 663 | fluxtop_dn(ig)=fluxtopvdn |
|---|
| 664 | |
|---|
| 665 | fluxtop_lw(ig) = real(nfluxtopi) |
|---|
| 666 | fluxabs_sw(ig) = real(-nfluxtopv) |
|---|
| 667 | fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD)) |
|---|
| 668 | fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD)) |
|---|
| 669 | |
|---|
| 670 | ! Flux absorbed by the surface. By MT2015. |
|---|
| 671 | fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig)) |
|---|
| 672 | |
|---|
| 673 | if(fluxtop_dn(ig).lt.0.0)then |
|---|
| 674 | print*,'Achtung! fluxtop_dn has lost the plot!' |
|---|
| 675 | print*,'fluxtop_dn=',fluxtop_dn(ig) |
|---|
| 676 | print*,'acosz=',acosz |
|---|
| 677 | print*,'temp= ',pt(ig,:) |
|---|
| 678 | print*,'pplay= ',pplay(ig,:) |
|---|
| 679 | call abort |
|---|
| 680 | endif |
|---|
| 681 | |
|---|
| 682 | ! Spectral output, for exoplanet observational comparison |
|---|
| 683 | if(specOLR)then |
|---|
| 684 | do nw=1,L_NSPECTI |
|---|
| 685 | OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth |
|---|
| 686 | end do |
|---|
| 687 | do nw=1,L_NSPECTV |
|---|
| 688 | !GSR_nu(ig,nw)=nfluxgndv_nu(nw) |
|---|
| 689 | OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth |
|---|
| 690 | end do |
|---|
| 691 | endif |
|---|
| 692 | |
|---|
| 693 | ! Finally, the heating rates |
|---|
| 694 | |
|---|
| 695 | DO l=2,L_NLAYRAD |
|---|
| 696 | dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) & |
|---|
| 697 | *gzlat(ig,L_NLAYRAD+1-l)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) |
|---|
| 698 | dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) & |
|---|
| 699 | *gzlat(ig,L_NLAYRAD+1-l)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) |
|---|
| 700 | END DO |
|---|
| 701 | |
|---|
| 702 | ! These are values at top of atmosphere |
|---|
| 703 | dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) & |
|---|
| 704 | *gzlat(ig,L_NLAYRAD)/(cpp*scalep*(plevrad(3)-plevrad(1))) |
|---|
| 705 | dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) & |
|---|
| 706 | *gzlat(ig,L_NLAYRAD)/(cpp*scalep*(plevrad(3)-plevrad(1))) |
|---|
| 707 | |
|---|
| 708 | |
|---|
| 709 | ! Optical thickness diagnostics (added by JVO) |
|---|
| 710 | if (diagdtau) then |
|---|
| 711 | do l=1,L_NLAYRAD |
|---|
| 712 | do nw=1,L_NSPECTV |
|---|
| 713 | int_dtauv(ig,l,nw) = 0.0d0 |
|---|
| 714 | DO k=1,L_NGAUSS |
|---|
| 715 | ! Output exp(-tau) because gweight ponderates exp and not tau itself |
|---|
| 716 | int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k) |
|---|
| 717 | ENDDO |
|---|
| 718 | enddo |
|---|
| 719 | do nw=1,L_NSPECTI |
|---|
| 720 | int_dtaui(ig,l,nw) = 0.0d0 |
|---|
| 721 | DO k=1,L_NGAUSS |
|---|
| 722 | ! Output exp(-tau) because gweight ponderates exp and not tau itself |
|---|
| 723 | int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k) |
|---|
| 724 | ENDDO |
|---|
| 725 | enddo |
|---|
| 726 | enddo |
|---|
| 727 | endif |
|---|
| 728 | |
|---|
| 729 | |
|---|
| 730 | ! Diagnostics : optical properties of haze and clouds (dtau, tau, k, w, g) : |
|---|
| 731 | do l = 2, L_LEVELS, 2 |
|---|
| 732 | lev2lay = L_NLAYRAD+1 - l/2 |
|---|
| 733 | ! Visible : |
|---|
| 734 | do nw = 1, L_NSPECTV |
|---|
| 735 | popthv(ig,lev2lay,nw,:) = 0.0d0 |
|---|
| 736 | popttv(ig,lev2lay,nw,:) = 0.0d0 |
|---|
| 737 | ! Optical thickness (dtau) : |
|---|
| 738 | popthv(ig,lev2lay,nw,1) = (diag_opthv(l,nw,1) + diag_opthv(l+1,nw,1)) / 2.0 |
|---|
| 739 | if (callclouds) then |
|---|
| 740 | popttv(ig,lev2lay,nw,1) = (diag_opttv(l,nw,1) + diag_opttv(l+1,nw,1)) / 2.0 |
|---|
| 741 | endif |
|---|
| 742 | ! Opacity (tau) : |
|---|
| 743 | do k = L_NLAYRAD, lev2lay, -1 |
|---|
| 744 | popthv(ig,lev2lay,nw,2) = popthv(ig,lev2lay,nw,2) + popthv(ig,k,nw,1) |
|---|
| 745 | if (callclouds) then |
|---|
| 746 | popttv(ig,lev2lay,nw,2) = popttv(ig,lev2lay,nw,2) + popttv(ig,k,nw,1) |
|---|
| 747 | endif |
|---|
| 748 | enddo |
|---|
| 749 | ! Extinction (k) : |
|---|
| 750 | popthv(ig,lev2lay,nw,3) = popthv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 751 | if (callclouds) then |
|---|
| 752 | popttv(ig,lev2lay,nw,3) = popttv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 753 | endif |
|---|
| 754 | ! Simple Scattering Albedo (w) and Asymmetry Parameter (g) : |
|---|
| 755 | popthv(ig,lev2lay,nw,4) = (diag_opthv(l,nw,2) + diag_opthv(l+1,nw,2)) / 2.0 |
|---|
| 756 | popthv(ig,lev2lay,nw,5) = (diag_opthv(l,nw,3) + diag_opthv(l+1,nw,3)) / 2.0 |
|---|
| 757 | if (callclouds) then |
|---|
| 758 | popttv(ig,lev2lay,nw,4) = (diag_opttv(l,nw,2) + diag_opttv(l+1,nw,2)) / 2.0 |
|---|
| 759 | popttv(ig,lev2lay,nw,5) = (diag_opttv(l,nw,3) + diag_opttv(l+1,nw,3)) / 2.0 |
|---|
| 760 | endif |
|---|
| 761 | ! DRAYAER, TAUGAS, DCONT : |
|---|
| 762 | popthv(ig,lev2lay,nw,6) = (diag_opthv(l,nw,4) + diag_opthv(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 763 | popthv(ig,lev2lay,nw,7) = (diag_opthv(l,nw,5) + diag_opthv(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 764 | popthv(ig,lev2lay,nw,8) = (diag_opthv(l,nw,6) + diag_opthv(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 765 | if (callclouds) then |
|---|
| 766 | popttv(ig,lev2lay,nw,6) = (diag_opttv(l,nw,4) + diag_opttv(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 767 | popttv(ig,lev2lay,nw,7) = (diag_opttv(l,nw,5) + diag_opttv(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 768 | popttv(ig,lev2lay,nw,8) = (diag_opttv(l,nw,6) + diag_opttv(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 769 | endif |
|---|
| 770 | enddo |
|---|
| 771 | ! Infra-Red |
|---|
| 772 | do nw=1,L_NSPECTI |
|---|
| 773 | popthi(ig,lev2lay,nw,:) = 0.0d0 |
|---|
| 774 | poptti(ig,lev2lay,nw,:) = 0.0d0 |
|---|
| 775 | ! Optical thickness (dtau) : |
|---|
| 776 | popthi(ig,lev2lay,nw,1) = (diag_opthi(l,nw,1) + diag_opthi(l+1,nw,1)) / 2.0 |
|---|
| 777 | if (callclouds) then |
|---|
| 778 | poptti(ig,lev2lay,nw,1) = (diag_optti(l,nw,1) + diag_optti(l+1,nw,1)) / 2.0 |
|---|
| 779 | endif |
|---|
| 780 | ! Opacity (tau) : |
|---|
| 781 | do k = L_NLAYRAD, lev2lay, -1 |
|---|
| 782 | popthi(ig,lev2lay,nw,2) = popthi(ig,lev2lay,nw,2) + popthi(ig,k,nw,1) |
|---|
| 783 | if (callclouds) then |
|---|
| 784 | poptti(ig,lev2lay,nw,2) = poptti(ig,lev2lay,nw,2) + poptti(ig,k,nw,1) |
|---|
| 785 | endif |
|---|
| 786 | enddo |
|---|
| 787 | ! Extinction (k) : |
|---|
| 788 | popthi(ig,lev2lay,nw,3) = popthi(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 789 | if (callclouds) then |
|---|
| 790 | poptti(ig,lev2lay,nw,3) = poptti(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 791 | endif |
|---|
| 792 | ! Simple Scattering Albedo (w) and Asymmetry Parameter (g) : |
|---|
| 793 | popthi(ig,lev2lay,nw,4) = (diag_opthi(l,nw,2) + diag_opthi(l+1,nw,2)) / 2.0 |
|---|
| 794 | popthi(ig,lev2lay,nw,5) = (diag_opthi(l,nw,3) + diag_opthi(l+1,nw,3)) / 2.0 |
|---|
| 795 | if (callclouds) then |
|---|
| 796 | poptti(ig,lev2lay,nw,4) = (diag_optti(l,nw,2) + diag_optti(l+1,nw,2)) / 2.0 |
|---|
| 797 | poptti(ig,lev2lay,nw,5) = (diag_optti(l,nw,3) + diag_optti(l+1,nw,3)) / 2.0 |
|---|
| 798 | endif |
|---|
| 799 | ! DRAYAER, TAUGAS, DCONT : |
|---|
| 800 | popthi(ig,lev2lay,nw,6) = (diag_opthi(l,nw,4) + diag_opthi(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 801 | popthi(ig,lev2lay,nw,7) = (diag_opthi(l,nw,5) + diag_opthi(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 802 | popthi(ig,lev2lay,nw,8) = (diag_opthi(l,nw,6) + diag_opthi(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 803 | if (callclouds) then |
|---|
| 804 | poptti(ig,lev2lay,nw,6) = (diag_optti(l,nw,4) + diag_optti(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 805 | poptti(ig,lev2lay,nw,7) = (diag_optti(l,nw,5) + diag_optti(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 806 | poptti(ig,lev2lay,nw,8) = (diag_optti(l,nw,6) + diag_optti(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay)) |
|---|
| 807 | endif |
|---|
| 808 | enddo |
|---|
| 809 | enddo |
|---|
| 810 | |
|---|
| 811 | |
|---|
| 812 | !----------------------------------------------------------------------- |
|---|
| 813 | end do ! End of big loop over every GCM column. |
|---|
| 814 | !----------------------------------------------------------------------- |
|---|
| 815 | |
|---|
| 816 | |
|---|
| 817 | !----------------------------------------------------------------------- |
|---|
| 818 | ! Additional diagnostics |
|---|
| 819 | !----------------------------------------------------------------------- |
|---|
| 820 | |
|---|
| 821 | ! IR spectral output, for exoplanet observational comparison |
|---|
| 822 | if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12 |
|---|
| 823 | |
|---|
| 824 | print*,'Saving scalar quantities in surf_vals.out...' |
|---|
| 825 | print*,'psurf = ', pplev(1,1),' Pa' |
|---|
| 826 | open(116,file='surf_vals.out') |
|---|
| 827 | write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), & |
|---|
| 828 | real(-nfluxtopv),real(nfluxtopi) |
|---|
| 829 | close(116) |
|---|
| 830 | |
|---|
| 831 | ! USEFUL COMMENT - Do Not Remove. |
|---|
| 832 | ! |
|---|
| 833 | ! if(specOLR)then |
|---|
| 834 | ! open(117,file='OLRnu.out') |
|---|
| 835 | ! do nw=1,L_NSPECTI |
|---|
| 836 | ! write(117,*) OLR_nu(1,nw) |
|---|
| 837 | ! enddo |
|---|
| 838 | ! close(117) |
|---|
| 839 | ! |
|---|
| 840 | ! open(127,file='OSRnu.out') |
|---|
| 841 | ! do nw=1,L_NSPECTV |
|---|
| 842 | ! write(127,*) OSR_nu(1,nw) |
|---|
| 843 | ! enddo |
|---|
| 844 | ! close(127) |
|---|
| 845 | ! endif |
|---|
| 846 | |
|---|
| 847 | ! OLR vs altitude: do it as a .txt file. |
|---|
| 848 | OLRz=.false. |
|---|
| 849 | if(OLRz)then |
|---|
| 850 | print*,'saving IR vertical flux for OLRz...' |
|---|
| 851 | open(118,file='OLRz_plevs.out') |
|---|
| 852 | open(119,file='OLRz.out') |
|---|
| 853 | do l=1,L_NLAYRAD |
|---|
| 854 | write(118,*) plevrad(2*l) |
|---|
| 855 | do nw=1,L_NSPECTI |
|---|
| 856 | write(119,*) fluxupi_nu(l,nw) |
|---|
| 857 | enddo |
|---|
| 858 | enddo |
|---|
| 859 | close(118) |
|---|
| 860 | close(119) |
|---|
| 861 | endif |
|---|
| 862 | |
|---|
| 863 | endif |
|---|
| 864 | |
|---|
| 865 | if (lastcall) then |
|---|
| 866 | IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi ) |
|---|
| 867 | IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv ) |
|---|
| 868 | IF( ALLOCATED( gasi_recomb ) ) DEALLOCATE( gasi_recomb ) |
|---|
| 869 | IF( ALLOCATED( gasv_recomb ) ) DEALLOCATE( gasv_recomb ) |
|---|
| 870 | IF( ALLOCATED( pqrold ) ) DEALLOCATE( pqrold ) |
|---|
| 871 | IF( ALLOCATED( useptold ) ) DEALLOCATE( useptold ) |
|---|
| 872 | !$OMP BARRIER |
|---|
| 873 | !$OMP MASTER |
|---|
| 874 | IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref ) |
|---|
| 875 | IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref ) |
|---|
| 876 | IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref ) |
|---|
| 877 | IF( ALLOCATED( gweight ) ) DEALLOCATE( gweight ) |
|---|
| 878 | !$OMP END MASTER |
|---|
| 879 | !$OMP BARRIER |
|---|
| 880 | endif |
|---|
| 881 | |
|---|
| 882 | |
|---|
| 883 | end subroutine callcorrk |
|---|
