| 1 | SUBROUTINE callradite(icount,ngrid,nlayer, |
|---|
| 2 | $ aerosol,albedo, |
|---|
| 3 | $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, |
|---|
| 4 | $ dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw) |
|---|
| 5 | |
|---|
| 6 | |
|---|
| 7 | IMPLICIT NONE |
|---|
| 8 | c======================================================================= |
|---|
| 9 | c subject: |
|---|
| 10 | c -------- |
|---|
| 11 | c Subroutine designed to call the main canonic |
|---|
| 12 | c radiative transfer subroutine "lwmain" et "swmain" |
|---|
| 13 | c to compute radiative heating and cooling rate and |
|---|
| 14 | c radiative fluxes to the surface. |
|---|
| 15 | c |
|---|
| 16 | c These calculations are only valid on the part of the atmosphere |
|---|
| 17 | c where Local Thermal Equilibrium (NLTE) is verified. In practice |
|---|
| 18 | c The calculations are only performed for the first "nlaylte" |
|---|
| 19 | c parameters (nlaylte is calculated by subroutine "nlthermeq" |
|---|
| 20 | c and stored in common "yomlw.h" |
|---|
| 21 | c |
|---|
| 22 | c |
|---|
| 23 | c |
|---|
| 24 | c The purpose of this subroutine is |
|---|
| 25 | c 1) Make some initial calculation at first call |
|---|
| 26 | c 2) Split the calculation in several sub-grid |
|---|
| 27 | c ("sub-domain") to save memory and |
|---|
| 28 | c be able run on a workstation at high resolution |
|---|
| 29 | c The sub-grid size is defined in dimradmars.h |
|---|
| 30 | c 3) call "lwmain" and "swmain" |
|---|
| 31 | c |
|---|
| 32 | c |
|---|
| 33 | c author: |
|---|
| 34 | c ------ |
|---|
| 35 | c Francois Forget / Christophe Hourdin |
|---|
| 36 | c |
|---|
| 37 | c This version modified to only calculate radiative tendencies over |
|---|
| 38 | c layers 1..NFLEV (set in dimradmars.h). Returns zero for higher |
|---|
| 39 | c layers, if any. |
|---|
| 40 | c In other routines, nlayermx -> nflev. |
|---|
| 41 | c Routines affected: lwflux, lwi, lwmain, lwxb, lwxd, lwxn. |
|---|
| 42 | c SRL 7/2000 |
|---|
| 43 | c |
|---|
| 44 | c definition: |
|---|
| 45 | c ---------- |
|---|
| 46 | c Here, solar band#1 is spectral interval between "long1vis" and "long2vis" |
|---|
| 47 | c set in dimradmars.h |
|---|
| 48 | c Here, solar band#2 is spectral interval between "long2vis" and "long3vis" |
|---|
| 49 | c set in dimradmars.h |
|---|
| 50 | c |
|---|
| 51 | c input: |
|---|
| 52 | c ----- |
|---|
| 53 | c icount counter of call to subroutine physic by gcm |
|---|
| 54 | c ngrid number of gridpoint of horizontal grid |
|---|
| 55 | c nlayer Number of layer |
|---|
| 56 | c |
|---|
| 57 | c aerosol(ngrid,nlayer,naerkind) aerosol extinction optical depth |
|---|
| 58 | c at reference wavelength "longrefvis" set |
|---|
| 59 | c in dimradmars.h , in each layer, for one of |
|---|
| 60 | c the "naerkind" kind of aerosol optical properties. |
|---|
| 61 | c |
|---|
| 62 | c albedo (ngrid,2) hemispheric surface albedo |
|---|
| 63 | c albedo (i,1) : mean albedo for solar band#1 |
|---|
| 64 | c (see below) |
|---|
| 65 | c albedo (i,2) : mean albedo for solar band#2 |
|---|
| 66 | c (see below) |
|---|
| 67 | c mu0(ngridmx) cos of solar zenith angle (=1 when sun at zenith) |
|---|
| 68 | c pplay(ngrid,nlayer) pressure (Pa) in the middle of each layer |
|---|
| 69 | c pplev(ngrid,nlayer+1) pressure (Pa) at boundaries of each layer |
|---|
| 70 | c pt(ngrid,nlayer) atmospheric temperature in each layer (K) |
|---|
| 71 | c tsurf(ngrid) surface temperature (K) |
|---|
| 72 | c fract(ngridmx) day fraction of the time interval |
|---|
| 73 | c =1 during the full day ; =0 during the night |
|---|
| 74 | c declin latitude of subsolar point |
|---|
| 75 | c dist_sol sun-Mars distance (AU) |
|---|
| 76 | c igout coordinate of analysed point for debugging |
|---|
| 77 | c |
|---|
| 78 | c output: |
|---|
| 79 | c ------- |
|---|
| 80 | c dtlw (ngrid,nlayer) longwave (IR) heating rate (K/s) |
|---|
| 81 | c dtsw(ngrid,nlayer) shortwave (Solar) heating rate (K/s) |
|---|
| 82 | c fluxsurf_lw(ngrid) surface downward flux tota LW (thermal IR) (W.m-2) |
|---|
| 83 | c fluxsurf_sw(ngrid,1) surface downward flux SW for solar band#1 (W.m-2) |
|---|
| 84 | c fluxsurf_sw(ngrid,2) surface downward flux SW for solar band#2 (W.m-2) |
|---|
| 85 | c |
|---|
| 86 | c fluxtop_lw(ngrid) outgoing upward flux tota LW (thermal IR) (W.m-2) |
|---|
| 87 | c fluxtop_sw(ngrid,1) outgoing upward flux SW for solar band#1 (W.m-2) |
|---|
| 88 | c fluxtop_sw(ngrid,2) outgoing upward flux SW for solar band#2 (W.m-2) |
|---|
| 89 | |
|---|
| 90 | c======================================================================= |
|---|
| 91 | c |
|---|
| 92 | c Declarations : |
|---|
| 93 | c ------------- |
|---|
| 94 | c |
|---|
| 95 | #include "dimensions.h" |
|---|
| 96 | #include "dimphys.h" |
|---|
| 97 | #include "dimradmars.h" |
|---|
| 98 | #include "comcstfi.h" |
|---|
| 99 | #include "callkeys.h" |
|---|
| 100 | #include "yomlw.h" |
|---|
| 101 | |
|---|
| 102 | |
|---|
| 103 | c----------------------------------------------------------------------- |
|---|
| 104 | c Input/Output |
|---|
| 105 | c ------------ |
|---|
| 106 | INTEGER Icount |
|---|
| 107 | INTEGER ngrid,nlayer |
|---|
| 108 | INTEGER igout |
|---|
| 109 | |
|---|
| 110 | REAL aerosol(ngrid,nlayer,naerkind) |
|---|
| 111 | REAL albedo(ngrid,2),emis(ngrid) |
|---|
| 112 | |
|---|
| 113 | REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer) |
|---|
| 114 | REAL pt(ngrid,nlayer) |
|---|
| 115 | REAL tsurf(ngrid) |
|---|
| 116 | REAL dist_sol,mu0(ngrid),fract(ngrid) |
|---|
| 117 | REAL dtlw(ngridmx,nlayermx),dtsw(ngridmx,nlayermx) |
|---|
| 118 | REAL fluxsurf_lw(ngridmx), fluxtop_lw(ngridmx) |
|---|
| 119 | REAL fluxsurf_sw(ngridmx,2), fluxtop_sw(ngridmx,2) |
|---|
| 120 | REAL flux(ngridmx,6) |
|---|
| 121 | |
|---|
| 122 | |
|---|
| 123 | c |
|---|
| 124 | c Local variables : |
|---|
| 125 | c ----------------- |
|---|
| 126 | |
|---|
| 127 | INTEGER l,ig, n |
|---|
| 128 | INTEGER jd,j,ig0,nd |
|---|
| 129 | |
|---|
| 130 | real cste_mars ! solar constant on Mars (Wm-2) |
|---|
| 131 | REAL ptlev(ngridmx,nlayermx+1) |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | INTEGER ndomain |
|---|
| 135 | parameter (ndomain = (ngridmx-1) / ndomainsz + 1) |
|---|
| 136 | |
|---|
| 137 | c Thermal IR net radiative budget (W m-2) |
|---|
| 138 | real znetrad(ndomainsz,nflev) |
|---|
| 139 | |
|---|
| 140 | real zfluxd_sw(ndomainsz,nflev+1,2) |
|---|
| 141 | real zfluxu_sw(ndomainsz,nflev+1,2) |
|---|
| 142 | |
|---|
| 143 | |
|---|
| 144 | REAL zplev(ndomainsz,nflev+1) |
|---|
| 145 | REAL zztlev(ndomainsz,nflev+1) |
|---|
| 146 | REAL zplay(ndomainsz,nflev) |
|---|
| 147 | REAL zt(ndomainsz,nflev) |
|---|
| 148 | REAL zaerosol(ndomainsz,nflev,naerkind) |
|---|
| 149 | REAL zalbedo(ndomainsz,2) |
|---|
| 150 | REAL zdp(ndomainsz,nflev) |
|---|
| 151 | REAL zdt0(ndomainsz) |
|---|
| 152 | |
|---|
| 153 | REAL zzdtlw(ndomainsz,nflev) |
|---|
| 154 | REAL zzdtsw(ndomainsz,nflev) |
|---|
| 155 | REAL zzflux(ndomainsz,6) |
|---|
| 156 | real zrmuz |
|---|
| 157 | |
|---|
| 158 | |
|---|
| 159 | c local saved variables |
|---|
| 160 | c --------------------- |
|---|
| 161 | |
|---|
| 162 | real pview(ngridmx) |
|---|
| 163 | save pview |
|---|
| 164 | |
|---|
| 165 | real zco2 ! volume fraction of CO2 in Mars atmosphere |
|---|
| 166 | DATA zco2/0.95/ |
|---|
| 167 | SAVE zco2 |
|---|
| 168 | |
|---|
| 169 | LOGICAL firstcall |
|---|
| 170 | DATA firstcall/.true./ |
|---|
| 171 | SAVE firstcall |
|---|
| 172 | |
|---|
| 173 | c---------------------------------------------------------------------- |
|---|
| 174 | |
|---|
| 175 | c Initialisation |
|---|
| 176 | c -------------- |
|---|
| 177 | |
|---|
| 178 | IF (firstcall) THEN |
|---|
| 179 | DO ig=1,ngrid |
|---|
| 180 | pview(ig)=1.66 ! cosecant of viewing angle |
|---|
| 181 | ENDDO |
|---|
| 182 | gcp = g/cpp |
|---|
| 183 | CALL SUAER |
|---|
| 184 | CALL SULW |
|---|
| 185 | |
|---|
| 186 | write(*,*) 'Splitting radiative calculations: ', |
|---|
| 187 | $ ' ngridmx,ngrid,ndomainsz,ndomain', |
|---|
| 188 | $ ngridmx,ngrid,ndomainsz,ndomain |
|---|
| 189 | if (ngridmx .EQ. 1) then |
|---|
| 190 | if (ndomainsz .NE. 1) then |
|---|
| 191 | print* |
|---|
| 192 | print*,'ATTENTION !!!' |
|---|
| 193 | print*,'pour tourner en 1D, ' |
|---|
| 194 | print*,'fixer ndomainsz=1 dans phymars/dimradmars.h' |
|---|
| 195 | print* |
|---|
| 196 | call exit(1) |
|---|
| 197 | endif |
|---|
| 198 | endif |
|---|
| 199 | firstcall=.false. |
|---|
| 200 | END IF |
|---|
| 201 | |
|---|
| 202 | c Starting loop on sub-domain |
|---|
| 203 | c ---------------------------- |
|---|
| 204 | |
|---|
| 205 | DO jd=1,ndomain |
|---|
| 206 | ig0=(jd-1)*ndomainsz |
|---|
| 207 | if (jd.eq.ndomain) then |
|---|
| 208 | nd=ngridmx-ig0 |
|---|
| 209 | else |
|---|
| 210 | nd=ndomainsz |
|---|
| 211 | endif |
|---|
| 212 | |
|---|
| 213 | |
|---|
| 214 | c Spliting input variable in sub-domain input variables |
|---|
| 215 | c --------------------------------------------------- |
|---|
| 216 | do l=1,nlaylte+1 |
|---|
| 217 | do ig = 1,nd |
|---|
| 218 | zplev(ig,l) = pplev(ig0+ig,l) |
|---|
| 219 | enddo |
|---|
| 220 | enddo |
|---|
| 221 | |
|---|
| 222 | do l=1,nlaylte |
|---|
| 223 | do ig = 1,nd |
|---|
| 224 | zplay(ig,l) = pplay(ig0+ig,l) |
|---|
| 225 | zt(ig,l) = pt(ig0+ig,l) |
|---|
| 226 | |
|---|
| 227 | c Thickness of each layer (Pa) : |
|---|
| 228 | zdp(ig,l)= pplev(ig0+ig,l) - pplev(ig0+ig,l+1) |
|---|
| 229 | enddo |
|---|
| 230 | enddo |
|---|
| 231 | |
|---|
| 232 | do n=1,naerkind |
|---|
| 233 | do l=1,nlaylte |
|---|
| 234 | do ig=1,nd |
|---|
| 235 | zaerosol(ig,l,n) = aerosol(ig0+ig,l,n) |
|---|
| 236 | enddo |
|---|
| 237 | enddo |
|---|
| 238 | enddo |
|---|
| 239 | |
|---|
| 240 | do j=1,2 |
|---|
| 241 | do ig = 1,nd |
|---|
| 242 | zalbedo(ig,j) = albedo(ig0+ig,j) |
|---|
| 243 | enddo |
|---|
| 244 | enddo |
|---|
| 245 | |
|---|
| 246 | c Intermediate levels: (computing tlev) |
|---|
| 247 | c --------------------------------------- |
|---|
| 248 | c Extrapolation for the air temperature above the surface |
|---|
| 249 | DO ig=1,nd |
|---|
| 250 | zztlev(ig,1)=zt(ig,1)+ |
|---|
| 251 | s (zplev(ig,1)-zplay(ig,1))* |
|---|
| 252 | s (zt(ig,1)-zt(ig,2))/(zplay(ig,1)-zplay(ig,2)) |
|---|
| 253 | |
|---|
| 254 | zdt0(ig) = tsurf(ig0+ig) - zztlev(ig,1) |
|---|
| 255 | ENDDO |
|---|
| 256 | |
|---|
| 257 | DO l=2,nlaylte |
|---|
| 258 | DO ig=1,nd |
|---|
| 259 | zztlev(ig,l)=0.5*(zt(ig,l-1)+zt(ig,l)) |
|---|
| 260 | ENDDO |
|---|
| 261 | ENDDO |
|---|
| 262 | |
|---|
| 263 | DO ig=1,nd |
|---|
| 264 | zztlev(ig,nlaylte+1)=zt(ig,nlaylte) |
|---|
| 265 | ENDDO |
|---|
| 266 | |
|---|
| 267 | |
|---|
| 268 | c Longwave ("lw") radiative transfer (= thermal infrared) |
|---|
| 269 | c ------------------------------------------------------- |
|---|
| 270 | |
|---|
| 271 | !PRINT*,'ig0',ig0 |
|---|
| 272 | !PRINT*,'icount',icount |
|---|
| 273 | !PRINT*,'nd',nd |
|---|
| 274 | !PRINT*,'nflev',nflev |
|---|
| 275 | !PRINT*,'zdp',zdp |
|---|
| 276 | !PRINT*,'zdt0',zdt0 |
|---|
| 277 | !PRINT*,'emis',emis(ig0+1) |
|---|
| 278 | !PRINT*,'zplev',zplev |
|---|
| 279 | !PRINT*,'zztlev',zztlev |
|---|
| 280 | !PRINT*,'zt',zt |
|---|
| 281 | !PRINT*,'zaerosol',zaerosol |
|---|
| 282 | !PRINT*,'zzdtlw',zzdtlw |
|---|
| 283 | !PRINT*,'fluxsurf_lw',fluxsurf_lw(ig0+1) |
|---|
| 284 | !PRINT*,'fluxtop_lw',fluxtop_lw(ig0+1) |
|---|
| 285 | !PRINT*,'znetrad',znetrad |
|---|
| 286 | |
|---|
| 287 | call lwmain (ig0,icount,nd,nflev |
|---|
| 288 | . ,zdp,zdt0,emis(ig0+1),zplev,zztlev,zt |
|---|
| 289 | . ,zaerosol,zzdtlw |
|---|
| 290 | . ,fluxsurf_lw(ig0+1),fluxtop_lw(ig0+1) |
|---|
| 291 | . ,znetrad) |
|---|
| 292 | |
|---|
| 293 | c Shortwave ("sw") radiative transfer (= solar radiation) |
|---|
| 294 | c ------------------------------------------------------- |
|---|
| 295 | c Mars solar constant (W m-2) |
|---|
| 296 | c 1370 W.m-2 is the solar constant at 1 AU. |
|---|
| 297 | cste_mars=1370./(dist_sol*dist_sol) |
|---|
| 298 | |
|---|
| 299 | call swmain ( nd, nflev, |
|---|
| 300 | S cste_mars, zalbedo, |
|---|
| 301 | S mu0(ig0+1), zdp, zplev, zaerosol, fract(ig0+1), |
|---|
| 302 | S zzdtsw, zfluxd_sw, zfluxu_sw) |
|---|
| 303 | |
|---|
| 304 | c Un-spliting output variable from sub-domain input variables |
|---|
| 305 | c ------------------------------------------------------------ |
|---|
| 306 | |
|---|
| 307 | do l=1,nlaylte |
|---|
| 308 | do ig = 1,nd |
|---|
| 309 | dtlw(ig0+ig,l) = zzdtlw(ig,l) |
|---|
| 310 | dtsw(ig0+ig,l) = zzdtsw(ig,l) |
|---|
| 311 | enddo |
|---|
| 312 | enddo |
|---|
| 313 | |
|---|
| 314 | do l=1,nlaylte+1 |
|---|
| 315 | do ig = 1,nd |
|---|
| 316 | ptlev(ig0+ig,l) = zztlev(ig,l) |
|---|
| 317 | enddo |
|---|
| 318 | enddo |
|---|
| 319 | |
|---|
| 320 | do ig = 1,nd |
|---|
| 321 | fluxsurf_sw(ig0+ig,1) = zfluxd_sw(ig,1,1) |
|---|
| 322 | fluxsurf_sw(ig0+ig,2) = zfluxd_sw(ig,1,2) |
|---|
| 323 | fluxtop_sw(ig0+ig,1) = zfluxu_sw(ig,nlaylte+1,1) |
|---|
| 324 | fluxtop_sw(ig0+ig,2) = zfluxu_sw(ig,nlaylte+1,2) |
|---|
| 325 | enddo |
|---|
| 326 | |
|---|
| 327 | |
|---|
| 328 | ENDDO ! (boucle jd=1, ndomain) |
|---|
| 329 | |
|---|
| 330 | |
|---|
| 331 | c Zero tendencies for any remaining layers between nlaylte and nlayer |
|---|
| 332 | if (nlayer.gt.nlaylte) then |
|---|
| 333 | do l = nlaylte+1, nlayer |
|---|
| 334 | do ig = 1, ngrid |
|---|
| 335 | dtlw(ig, l) = 0. |
|---|
| 336 | dtsw(ig, l) = 0. |
|---|
| 337 | enddo |
|---|
| 338 | enddo |
|---|
| 339 | endif |
|---|
| 340 | |
|---|
| 341 | |
|---|
| 342 | c Output for debugging if lwrite=T |
|---|
| 343 | c -------------------------------- |
|---|
| 344 | c Write all nlayer layers, even though only nlaylte layers may have |
|---|
| 345 | c non-zero tendencies. |
|---|
| 346 | |
|---|
| 347 | IF(lwrite) THEN |
|---|
| 348 | PRINT*,'Diagnotique for the radiation' |
|---|
| 349 | PRINT*,'albedo, emissiv, mu0,fract,fluxsurf_lw,fluxsurf_sw' |
|---|
| 350 | PRINT*,albedo(igout,1),emis(igout),mu0(igout), |
|---|
| 351 | s fract(igout), fluxsurf_lw(igout), |
|---|
| 352 | $ fluxsurf_sw(igout,1)+fluxsurf_sw(igout,2) |
|---|
| 353 | PRINT*,'Tlay Tlev Play Plev dT/dt SW dT/dt LW (K/s)' |
|---|
| 354 | PRINT*,'daysec',daysec |
|---|
| 355 | DO l=1,nlayer |
|---|
| 356 | PRINT*,pt(igout,l),ptlev(igout,l), |
|---|
| 357 | s pplay(igout,l),pplev(igout,l), |
|---|
| 358 | s dtsw(igout,l),dtlw(igout,l) |
|---|
| 359 | ENDDO |
|---|
| 360 | ENDIF |
|---|
| 361 | |
|---|
| 362 | |
|---|
| 363 | return |
|---|
| 364 | end |
|---|
| 365 | |
|---|
