Changeset 1483
- Timestamp:
- Oct 14, 2015, 8:56:32 PM (9 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 2 edited
-
README (modified) (1 diff)
-
libf/phystd/callcorrk.F90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r1482 r1483 1091 1091 - CO2 Ice Albedo North/South dichotomy has been removed. CO2 ice Albedo is no longer stocked in start file. You can now 1092 1092 prescribe CO2 ice albedo in callphys.def with "albedoco2ice". Its default value is 0.5. 1093 1094 == 14/10/2015 == MT 1095 - Cleanup of callcorrk.F90 -
trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
r1482 r1483 13 13 use watercommon_h 14 14 use datafile_mod, only: datadir 15 ! use ioipsl_getincom16 15 use ioipsl_getincom_p 17 16 use gases_h … … 21 20 use comcstfi_mod, only: pi, mugaz, cpp 22 21 use callkeys_mod, only: varactive,diurnal,tracer,water,nosurf,varfixed,satval, & 23 kastprof,strictboundcorrk,specOLR,CLFvarying22 kastprof,strictboundcorrk,specOLR,CLFvarying 24 23 25 24 implicit none … … 44 43 ! Layer #1 is the layer near the ground. 45 44 ! Layer #nlayer is the layer at the top. 46 47 INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns 48 INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers 49 REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air) 50 integer,intent(in) :: nq ! number of tracers 51 REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg.m-2) 52 REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral SW albedo.MT2015. 53 REAL,INTENT(IN) :: emis(ngrid) ! LW emissivity 54 real,intent(in) :: mu0(ngrid) ! cosine of sun incident angle 55 REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) 56 REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa) 57 REAL,INTENT(IN) :: pt(ngrid,nlayer) ! air temperature (K) 58 REAL,INTENT(IN) :: tsurf(ngrid) ! surface temperature (K) 59 REAL,INTENT(IN) :: fract(ngrid) ! fraction of day 60 REAL,INTENT(IN) :: dist_star ! distance star-planet (AU) 61 REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol tau (kg/kg) 62 real,intent(in) :: muvar(ngrid,nlayer+1) 63 REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! heating rate (K/s) due to LW 64 REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! heating rate (K/s) due to SW 65 REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! incident LW flux to surf (W/m2) 66 REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! incident SW flux to surf (W/m2) 67 REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! absorbed SW flux by the surface (W/m2). By MT2015. 68 REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! outgoing LW flux to space (W/m2) 69 REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by planet (W/m2) 70 REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! incident top of atmosphere SW flux (W/m2) 71 REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1) 72 REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1) 73 REAL,INTENT(OUT) :: tau_col(ngrid) ! diagnostic from aeropacity 74 REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015 75 ! for H2O cloud fraction in aeropacity 76 real,intent(in) :: cloudfrac(ngrid,nlayer) 77 real,intent(out) :: totcloudfrac(ngrid) 45 !----------------------------------------------------------------------- 46 47 48 ! INPUT 49 INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns. 50 INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers. 51 REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air). 52 INTEGER,INTENT(IN) :: nq ! Number of tracers. 53 REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2). 54 REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015 55 REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity. 56 REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle. 57 REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa). 58 REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa). 59 REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K). 60 REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K). 61 REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day. 62 REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU). 63 REAL,INTENT(IN) :: muvar(ngrid,nlayer+1) 64 REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer) ! Fraction of clouds (%). 78 65 logical,intent(in) :: clearsky 79 logical,intent(in) :: firstcall ! signals first call to physics 80 logical,intent(in) :: lastcall ! signals last call to physics 81 82 ! Globally varying aerosol optical properties on GCM grid 83 ! Not needed everywhere so not in radcommon_h 66 logical,intent(in) :: firstcall ! Signals first call to physics. 67 logical,intent(in) :: lastcall ! Signals last call to physics. 68 69 ! OUTPUT 70 REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg). 71 REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation. 72 REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation. 73 REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2). 74 REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2) 75 REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015. 76 REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2). 77 REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2). 78 REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2). 79 REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1). 80 REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1). 81 REAL,INTENT(OUT) :: tau_col(ngrid) ! Diagnostic from aeropacity. 82 REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015 83 REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! Column Fraction of clouds (%). 84 85 86 87 88 89 ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h. 84 90 REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind) 85 91 REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind) 86 92 REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind) 87 88 93 REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind) 89 94 REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind) … … 93 98 ! REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these... 94 99 95 REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)100 REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m) 96 101 REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance 97 102 !$OMP THREADPRIVATE(reffrad,nueffrad) … … 99 104 !----------------------------------------------------------------------- 100 105 ! Declaration of the variables required by correlated-k subroutines 101 ! Numbered from top to bottom unlike in the GCM! 106 ! Numbered from top to bottom (unlike in the GCM) 107 !----------------------------------------------------------------------- 102 108 103 109 REAL*8 tmid(L_LEVELS),pmid(L_LEVELS) 104 110 REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS) 105 111 106 !Optical values for the optci/cv subroutines112 ! Optical values for the optci/cv subroutines 107 113 REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV) 108 114 REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS) … … 118 124 REAL*8 tauaero(L_LEVELS+1,naerkind) 119 125 REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn 120 real*8 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2)121 real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2)122 real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic126 REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2). 127 REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2). 128 REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic. 123 129 REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD) 124 130 REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD) 125 131 REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD) 126 132 REAL*8 albi,acosz 127 REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo. 128 129 INTEGER ig,l,k,nw,iaer,irad 130 INTEGER icount 133 REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo. 134 135 INTEGER ig,l,k,nw,iaer 131 136 132 137 real szangle … … 136 141 real*8 taugsurf(L_NSPECTV,L_NGAUSS-1) 137 142 real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1) 138 139 real*8 qvar(L_LEVELS) ! mixing ratio of variable component (mol/mol) 140 141 ! Local aerosol optical properties for each column on RADIATIVE grid 143 real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol). 144 145 ! Local aerosol optical properties for each column on RADIATIVE grid. 142 146 real*8,save :: QXVAER(L_LEVELS+1,L_NSPECTV,naerkind) 143 147 real*8,save :: QSVAER(L_LEVELS+1,L_NSPECTV,naerkind) … … 155 159 156 160 157 ! Misc. 158 logical nantest 159 real*8 tempv(L_NSPECTV) 160 real*8 tempi(L_NSPECTI) 161 ! Miscellaneous : 161 162 real*8 temp,temp1,temp2,pweight 162 163 character(len=10) :: tmp1 163 164 character(len=10) :: tmp2 164 165 165 ! for fixed water vapour profiles 166 ! For fixed water vapour profiles. 166 167 integer i_var 167 168 real RH 168 169 real*8 pq_temp(nlayer) 170 ! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!! 169 171 real ptemp, Ttemp, qsat 170 172 171 ! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!172 173 !real ptime, pday174 173 logical OLRz 175 174 real*8 NFLUXGNDV_nu(L_NSPECTV) 176 175 177 178 ! for weird cloud test 179 real pqtest(ngrid,nlayer,nq) 180 181 real maxrad, minrad 182 183 real,external :: CBRT 184 185 ! included by RW for runaway greenhouse 1D study 176 ! Included by RW for runaway greenhouse 1D study. 186 177 real vtmp(nlayer) 187 178 REAL*8 muvarrad(L_LEVELS) 188 179 189 ! included by MT for albedo calculations.180 ! Included by MT for albedo calculations. 190 181 REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation. 191 182 183 184 192 185 !=============================================================== 193 ! Initialization on first call 186 ! I.a Initialization on first call 187 !=============================================================== 188 194 189 195 190 qxvaer(:,:,:)=0.0 … … 203 198 if(firstcall) then 204 199 205 !!! ALLOCATED instances are necessary because of CLFvarying 206 !!! strategy to call callcorrk twice in physiq... 200 !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...) 207 201 IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind)) 208 202 IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind)) … … 219 213 call su_aer_radii(ngrid,nlayer,reffrad,nueffrad) 220 214 215 221 216 !-------------------------------------------------- 222 ! set up correlated k 217 ! Set up correlated k 218 !-------------------------------------------------- 219 220 223 221 print*, "callcorrk: Correlated-k data base folder:",trim(datadir) 224 222 call getin_p("corrkdir",corrkdir) … … 229 227 banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)) 230 228 231 call setspi ! basic infrared properties232 call setspv ! basic visible properties233 call sugas_corrk ! set up gaseous absorption properties234 call suaer_corrk ! set up aerosol optical properties229 call setspi ! Basic infrared properties. 230 call setspv ! Basic visible properties. 231 call sugas_corrk ! Set up gaseous absorption properties. 232 call suaer_corrk ! Set up aerosol optical properties. 235 233 236 234 … … 244 242 245 243 if (ngrid.eq.1) then 246 PRINT*, 'Simulate global averaged conditions ?'247 global1d = .false. ! default value248 call getin_p("global1d",global1d)249 write(*,*) "global1d = ",global1d250 ! Test of incompatibility:251 ! if global1d is true, there should not be any diurnal cycle252 if (global1d.and.diurnal) then253 print*,'if global1d is true, diurnal must be set to false'254 stop255 endif256 257 if (global1d) then258 PRINT *,'Solar Zenith angle (deg.) ?'259 PRINT *,'(assumed for averaged solar flux S/4)'260 szangle=60.0 ! default value261 call getin_p("szangle",szangle)262 write(*,*) "szangle = ",szangle263 endif244 PRINT*, 'Simulate global averaged conditions ?' 245 global1d = .false. ! default value 246 call getin_p("global1d",global1d) 247 write(*,*) "global1d = ",global1d 248 249 ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle. 250 if (global1d.and.diurnal) then 251 print*,'if global1d is true, diurnal must be set to false' 252 stop 253 endif 254 255 if (global1d) then 256 PRINT *,'Solar Zenith angle (deg.) ?' 257 PRINT *,'(assumed for averaged solar flux S/4)' 258 szangle=60.0 ! default value 259 call getin_p("szangle",szangle) 260 write(*,*) "szangle = ",szangle 261 endif 264 262 endif 265 263 266 264 end if ! of if (firstcall) 267 265 266 268 267 !======================================================================= 269 ! Initialization on every call 270 268 ! I.b Initialization on every call 269 !======================================================================= 270 271 271 272 !-------------------------------------------------- 272 273 ! Effective radius and variance of the aerosols 274 !-------------------------------------------------- 275 273 276 do iaer=1,naerkind 274 277 275 if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.278 if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! Treat condensed co2 particles. 276 279 call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2)) 277 280 print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um' 278 281 print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um' 279 282 end if 280 if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols 283 284 if ((iaer.eq.iaero_h2o).and.water) then ! Treat condensed water particles. To be generalized for other aerosols ... 281 285 call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, & 282 286 reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o)) … … 284 288 print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um' 285 289 endif 290 286 291 if(iaer.eq.iaero_dust)then 287 292 call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust)) 288 293 print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um' 289 294 endif 295 290 296 if(iaer.eq.iaero_h2so4)then 291 297 call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4)) 292 298 print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um' 293 299 endif 300 294 301 if(iaer.eq.iaero_back2lay)then 295 302 call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev) 296 303 endif 297 end do !iaer=1,naerkind 298 299 300 ! how much light we get 304 305 end do !iaer=1,naerkind. 306 307 308 ! How much light do we get ? 301 309 do nw=1,L_NSPECTV 302 310 stel(nw)=stellarf(nw)/(dist_star**2) 303 311 end do 304 312 313 ! Get 3D aerosol optical properties. 305 314 call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, & 306 315 QVISsQREF3d,omegaVIS3d,gVIS3d, & 307 316 QIRsQREF3d,omegaIR3d,gIR3d, & 308 QREFvis3d,QREFir3d) ! get 3D aerosol optical properties 309 317 QREFvis3d,QREFir3d) 318 319 ! Get aerosol optical depths. 310 320 call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, & 311 321 reffrad,QREFvis3d,QREFir3d, & 312 tau_col,cloudfrac,totcloudfrac,clearsky) ! get aerosol optical depths322 tau_col,cloudfrac,totcloudfrac,clearsky) 313 323 314 !----------------------------------------------------------------------- 315 ! Starting Big Loop over every GCM column 316 do ig=1,ngrid 324 325 326 !----------------------------------------------------------------------- 327 do ig=1,ngrid ! Starting Big Loop over every GCM column 328 !----------------------------------------------------------------------- 329 317 330 318 331 !======================================================================= 319 ! Transformation of the GCM variables 320 321 !----------------------------------------------------------------------- 322 ! Aerosol optical properties Qext, Qscat and g 323 ! The transformation in the vertical is the same as for temperature 332 ! II. Transformation of the GCM variables 333 !======================================================================= 334 335 336 !----------------------------------------------------------------------- 337 ! Aerosol optical properties Qext, Qscat and g. 338 ! The transformation in the vertical is the same as for temperature. 339 !----------------------------------------------------------------------- 324 340 325 ! shortwave 341 326 342 do iaer=1,naerkind 327 DO nw=1,L_NSPECTV 343 ! Shortwave. 344 do nw=1,L_NSPECTV 345 328 346 do l=1,nlayer 329 347 … … 349 367 gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 350 368 351 end do 369 end do ! nlayer 352 370 353 371 qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer) … … 360 378 gvaer(2*nlayer+1,nw,iaer)=0. 361 379 362 end do 363 364 ! longwave 365 DO nw=1,L_NSPECTI380 end do ! L_NSPECTV 381 382 do nw=1,L_NSPECTI 383 ! Longwave 366 384 do l=1,nlayer 367 385 … … 387 405 giaer(2*l+1,nw,iaer)=(temp1+temp2)/2 388 406 389 end do 407 end do ! nlayer 390 408 391 409 qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer) … … 398 416 giaer(2*nlayer+1,nw,iaer)=0. 399 417 400 end do 401 end do 402 403 ! test / correct for freaky s. s. albedo values 418 end do ! L_NSPECTI 419 420 end do ! naerkind 421 422 ! Test / Correct for freaky s. s. albedo values. 404 423 do iaer=1,naerkind 405 424 do k=1,L_LEVELS+1 … … 427 446 end do 428 447 429 end do 430 end do 448 end do ! L_LEVELS 449 end do ! naerkind 431 450 432 451 !----------------------------------------------------------------------- 433 452 ! Aerosol optical depths 453 !----------------------------------------------------------------------- 434 454 435 455 do iaer=1,naerkind ! a bug was here … … 437 457 438 458 pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ & 439 (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1)) 440 459 (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1)) 441 460 temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer) 442 443 461 tauaero(2*k+2,iaer)=max(temp*pweight,0.d0) 444 462 tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0) 445 ! 463 446 464 end do 447 465 ! boundary conditions … … 450 468 !tauaero(1,iaer) = 0. 451 469 !tauaero(L_LEVELS+1,iaer) = 0. 452 end do 453 454 ! Albedo and emissivity 455 albi=1-emis(ig) ! Longwave. 456 DO nw=1,L_NSPECTV ! Shortwave loop. 470 471 end do ! naerkind 472 473 ! Albedo and Emissivity. 474 albi=1-emis(ig) ! Long Wave. 475 DO nw=1,L_NSPECTV ! Short Wave loop. 457 476 albv(nw)=albedo(ig,nw) 458 477 ENDDO 459 478 460 if (nosurf) then 479 if (nosurf) then ! Case with no surface. 461 480 DO nw=1,L_NSPECTV 462 481 if(albv(nw).gt.0.0) then … … 468 487 endif 469 488 470 if ((ngrid.eq.1).and.(global1d)) then ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight489 if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight. 471 490 acosz = cos(pi*szangle/180.0) 472 491 print*,'acosz=',acosz,', szangle=',szangle 473 492 else 474 acosz=mu0(ig) ! cosine of sun incident angle : 3D simulations or local 1D simulations using latitude493 acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude. 475 494 endif 476 495 477 !!! JL13: in the following, I changed some indices in the interpolations so that the model rsults are less dependent on the number of layers 478 !!! the older verions are commented with the commetn !JL13index 479 480 481 !----------------------------------------------------------------------- 482 ! Water vapour (to be generalised for other gases eventually) 496 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 497 !!! Note by JL13 : In the following, some indices were changed in the interpolations, 498 !!! so that the model results are less dependent on the number of layers ! 499 !!! 500 !!! --- The older versions are commented with the comment !JL13index --- 501 !!! 502 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 503 504 505 !----------------------------------------------------------------------- 506 ! Water vapour (to be generalised for other gases eventually ...) 507 !----------------------------------------------------------------------- 483 508 484 509 if(varactive)then … … 495 520 elseif(varfixed)then 496 521 497 do l=1,nlayer ! here we will assign fixed water vapour profiles globally522 do l=1,nlayer ! Here we will assign fixed water vapour profiles globally. 498 523 RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98) 499 524 if(RH.lt.0.0) RH=0.0 … … 511 536 qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2 512 537 end do 538 513 539 qvar(1)=qvar(2) 514 540 … … 527 553 qvar(k) = 1.0D-7 528 554 end do 529 end if 555 end if ! varactive/varfixed 530 556 531 557 if(.not.kastprof)then 532 ! IMPORTANT: Now convert from kg/kg to mol/mol558 ! IMPORTANT: Now convert from kg/kg to mol/mol. 533 559 do k=1,L_LEVELS 534 560 qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi)) … … 537 563 538 564 !----------------------------------------------------------------------- 539 ! kcm mode only 565 ! kcm mode only ! 566 !----------------------------------------------------------------------- 567 540 568 if(kastprof)then 541 569 542 ! initial values equivalent to mugaz570 ! Initial values equivalent to mugaz. 543 571 DO l=1,nlayer 544 572 muvarrad(2*l) = mugaz … … 549 577 550 578 DO l=1,nlayer 551 muvarrad(2*l) = muvar(ig,nlayer+2-l)579 muvarrad(2*l) = muvar(ig,nlayer+2-l) 552 580 muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + & 553 muvar(ig,max(nlayer+1-l,1)))/2581 muvar(ig,max(nlayer+1-l,1)))/2 554 582 END DO 555 583 556 584 muvarrad(1) = muvarrad(2) 557 muvarrad(2*nlayer+1) =muvar(ig,1)585 muvarrad(2*nlayer+1) = muvar(ig,1) 558 586 559 587 print*,'Recalculating qvar with VARIABLE epsi for kastprof' 560 588 print*,'Assumes that the variable gas is H2O!!!' 561 589 print*,'Assumes that there is only one tracer' 590 562 591 !i_var=igcm_h2o_vap 563 592 i_var=1 593 564 594 if(nq.gt.1)then 565 595 print*,'Need 1 tracer only to run kcm1d.e' 566 596 stop 567 597 endif 598 568 599 do l=1,nlayer 569 600 vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi)) … … 580 611 print*,'Warning: reducing qvar in callcorrk.' 581 612 print*,'Temperature profile no longer consistent ', & 582 'with saturated H2O. qsat=',satval 613 'with saturated H2O. qsat=',satval 614 583 615 do k=1,L_LEVELS 584 616 qvar(k) = qvar(k)*satval … … 594 626 muvarrad(1) = muvarrad(2) 595 627 muvarrad(2*nlayer+1)=muvar(ig,1) 596 endif 597 598 ! Keep values inside limits for which we have radiative transfer coefficients 599 if(L_REFVAR.gt.1)then ! there was a bug here!628 endif ! if kastprof 629 630 ! Keep values inside limits for which we have radiative transfer coefficients !!! 631 if(L_REFVAR.gt.1)then ! (there was a bug here) 600 632 do k=1,L_LEVELS 601 633 if(qvar(k).lt.wrefvar(1))then … … 609 641 !----------------------------------------------------------------------- 610 642 ! Pressure and temperature 643 !----------------------------------------------------------------------- 611 644 612 645 DO l=1,nlayer … … 618 651 619 652 plevrad(1) = 0. 620 plevrad(2) = 0. !! Trick to have correct calculations of fluxes in gflux(i/v).F, 621 !! but the pmid levels are not impacted by this change. 653 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. 622 654 623 655 tlevrad(1) = tlevrad(2) … … 653 685 ! tmid(L_LEVELS) = tlevrad(L_LEVELS-1) 654 686 655 ! test for out-of-bounds pressure687 ! Test for out-of-bounds pressure. 656 688 if(plevrad(3).lt.pgasmin)then 657 689 print*,'Minimum pressure is outside the radiative' … … 664 696 endif 665 697 666 ! test for out-of-bounds temperature698 ! Test for out-of-bounds temperature. 667 699 do k=1,L_LEVELS 668 700 if(tlevrad(k).lt.tgasmin)then … … 733 765 734 766 !======================================================================= 735 ! Calling the main radiative transfer subroutines 736 737 738 Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed 767 ! III. Calling the main radiative transfer subroutines 768 !======================================================================= 769 770 771 Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed. 739 772 glat_ig=glat(ig) 740 773 741 774 !----------------------------------------------------------------------- 742 ! Shortwave 743 744 if(fract(ig) .ge. 1.0e-4) then ! only during daylight! 775 ! Short Wave Part 776 !----------------------------------------------------------------------- 777 778 if(fract(ig) .ge. 1.0e-4) then ! Only during daylight. 745 779 if((ngrid.eq.1).and.(global1d))then 746 780 do nw=1,L_NSPECTV 747 stel_fract(nw)= stel(nw)* 0.25 / acosz 748 ! globally averaged = divide by 4 749 ! but we correct for solar zenith angle 781 stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle 750 782 end do 751 783 else … … 753 785 stel_fract(nw)= stel(nw) * fract(ig) 754 786 end do 755 endif 787 endif 788 756 789 call optcv(dtauv,tauv,taucumv,plevrad, & 757 790 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, & … … 763 796 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf) 764 797 765 else ! during the night, fluxes = 0798 else ! During the night, fluxes = 0. 766 799 nfluxtopv = 0.0d0 767 800 fluxtopvdn = 0.0d0 … … 789 822 790 823 !----------------------------------------------------------------------- 791 ! Longwave 824 ! Long Wave Part 825 !----------------------------------------------------------------------- 792 826 793 827 call optci(plevrad,tlevrad,dtaui,taucumi, & … … 850 884 *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1))) 851 885 852 ! --------------------------------------------------------------- 853 end do ! end of big loop over every GCM column (ig = 1:ngrid) 886 887 !----------------------------------------------------------------------- 888 end do ! End of big loop over every GCM column. 889 !----------------------------------------------------------------------- 890 854 891 855 892 856 893 !----------------------------------------------------------------------- 857 894 ! Additional diagnostics 858 859 ! IR spectral output, for exoplanet observational comparison 860 861 895 !----------------------------------------------------------------------- 896 897 ! IR spectral output, for exoplanet observational comparison 862 898 if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12 863 899 864 print*,'Saving scalar quantities in surf_vals.out...' 865 print*,'psurf = ', pplev(1,1),' Pa' 866 open(116,file='surf_vals.out') 867 write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), & 868 real(-nfluxtopv),real(nfluxtopi) 869 close(116) 870 871 ! I am useful, please don`t remove me! 900 print*,'Saving scalar quantities in surf_vals.out...' 901 print*,'psurf = ', pplev(1,1),' Pa' 902 open(116,file='surf_vals.out') 903 write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), & 904 real(-nfluxtopv),real(nfluxtopi) 905 close(116) 906 907 908 ! USEFUL COMMENT - Do Not Remove. 909 ! 872 910 ! if(specOLR)then 873 911 ! open(117,file='OLRnu.out') … … 884 922 ! endif 885 923 886 ! OLR vs altitude: do it as a .txt file 887 888 889 890 891 892 893 894 895 896 897 898 899 900 924 ! OLR vs altitude: do it as a .txt file. 925 OLRz=.false. 926 if(OLRz)then 927 print*,'saving IR vertical flux for OLRz...' 928 open(118,file='OLRz_plevs.out') 929 open(119,file='OLRz.out') 930 do l=1,L_NLAYRAD 931 write(118,*) plevrad(2*l) 932 do nw=1,L_NSPECTI 933 write(119,*) fluxupi_nu(l,nw) 934 enddo 935 enddo 936 close(118) 937 close(119) 938 endif 901 939 902 940 endif 903 941 904 ! see physiq.F for explanations about CLFvarying. This is temporary.942 ! See physiq.F for explanations about CLFvarying. This is temporary. 905 943 if (lastcall .and. .not.CLFvarying) then 906 944 IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
Note: See TracChangeset
for help on using the changeset viewer.
