Changeset 1482 for trunk/LMDZ.GENERIC/libf/phystd
- Timestamp:
- Oct 14, 2015, 5:05:47 PM (9 years ago)
- Location:
- trunk/LMDZ.GENERIC/libf/phystd
- Files:
-
- 11 edited
-
callcorrk.F90 (modified) (10 diffs)
-
callkeys_mod.F90 (modified) (1 diff)
-
condense_co2.F90 (modified) (6 diffs)
-
hydrol.F90 (modified) (10 diffs)
-
inifis.F (modified) (1 diff)
-
phyredem.F90 (modified) (2 diffs)
-
physiq.F90 (modified) (18 diffs)
-
sfluxv.F (modified) (5 diffs)
-
surfdat_h.F90 (modified) (2 diffs)
-
surfini.F (modified) (1 diff)
-
tabfi.F (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
r1423 r1482 1 1 subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, & 2 albedo, emis,mu0,pplev,pplay,pt,&2 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 3 3 tsurf,fract,dist_star,aerosol,muvar, & 4 4 dtlw,dtsw,fluxsurf_lw, & 5 fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn, & 5 fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, & 6 fluxabs_sw,fluxtop_dn, & 6 7 OLR_nu,OSR_nu, & 7 8 tau_col,cloudfrac,totcloudfrac, & … … 49 50 integer,intent(in) :: nq ! number of tracers 50 51 REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg.m-2) 51 REAL,INTENT(IN) :: albedo(ngrid ) ! SW albedo52 REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral SW albedo.MT2015. 52 53 REAL,INTENT(IN) :: emis(ngrid) ! LW emissivity 53 54 real,intent(in) :: mu0(ngrid) ! cosine of sun incident angle … … 64 65 REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! incident LW flux to surf (W/m2) 65 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. 66 68 REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! outgoing LW flux to space (W/m2) 67 69 REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by planet (W/m2) … … 70 72 REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1) 71 73 REAL,INTENT(OUT) :: tau_col(ngrid) ! diagnostic from aeropacity 74 REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015 72 75 ! for H2O cloud fraction in aeropacity 73 76 real,intent(in) :: cloudfrac(ngrid,nlayer) … … 121 124 REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD) 122 125 REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD) 123 REAL*8 albi,albv,acosz 126 REAL*8 albi,acosz 127 REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo. 124 128 125 129 INTEGER ig,l,k,nw,iaer,irad … … 182 186 real vtmp(nlayer) 183 187 REAL*8 muvarrad(L_LEVELS) 188 189 ! included by MT for albedo calculations. 190 REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation. 184 191 185 192 !=============================================================== … … 212 219 call su_aer_radii(ngrid,nlayer,reffrad,nueffrad) 213 220 214 215 216 221 !-------------------------------------------------- 217 222 ! set up correlated k … … 448 453 449 454 ! Albedo and emissivity 450 albi=1-emis(ig) ! longwave 451 albv=albedo(ig) ! shortwave 452 453 if(nosurf.and.(albv.gt.0.0))then 454 print*,'For open lower boundary in callcorrk must' 455 print*,'have surface albedo set to zero!' 456 call abort 457 endif 455 albi=1-emis(ig) ! Longwave. 456 DO nw=1,L_NSPECTV ! Shortwave loop. 457 albv(nw)=albedo(ig,nw) 458 ENDDO 459 460 if (nosurf) then 461 DO nw=1,L_NSPECTV 462 if(albv(nw).gt.0.0) then 463 print*,'For open lower boundary in callcorrk must' 464 print*,'have spectral surface band albedos all set to zero!' 465 call abort 466 endif 467 ENDDO 468 endif 458 469 459 470 if ((ngrid.eq.1).and.(global1d)) then ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight … … 764 775 end if 765 776 777 778 ! Equivalent Albedo Calculation. MT2015 779 if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight. 780 DO nw=1,L_NSPECTV 781 albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw) 782 ENDDO 783 albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/sum(nfluxgndv_nu(1:L_NSPECTV)) 784 albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV)) 785 else 786 albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0. 787 endif 788 789 766 790 !----------------------------------------------------------------------- 767 791 ! Longwave … … 786 810 fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD)) 787 811 fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD)) 812 813 ! Flux absorbed by the surface. By MT2015. 814 fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig)) 788 815 789 816 if(fluxtop_dn(ig).lt.0.0)then -
trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
r1423 r1482 66 66 real pceil 67 67 real albedosnow 68 real albedoco2ice 68 69 real maxicethick 69 70 real Tsaldiff -
trunk/LMDZ.GENERIC/libf/phystd/condense_co2.F90
r1477 r1482 1 1 subroutine condense_co2(ngrid,nlayer,nq,ptimestep, & 2 pcapcal,pplay,pplev,ptsrf,pt, & 3 pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, & 4 piceco2,psolaralb,pemisurf, & 5 pdtc,pdtsrfc,pdpsrf,pduc,pdvc, & 2 pcapcal,pplay,pplev,ptsrf,pt, & 3 pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, & 4 piceco2,albedo,pemisurf, & 5 albedo_bareground,albedo_co2_ice_SPECTV, & 6 pdtc,pdtsrfc,pdpsrf,pduc,pdvc, & 6 7 pdqc) 7 8 8 use radinc_h, only : naerkind9 use radinc_h, only : L_NSPECTV, naerkind 9 10 use gases_h, only: gfrac, igas_co2 10 11 use radii_mod, only : co2_reffrad 11 12 use aerosol_mod, only : iaero_co2 12 USE surfdat_h, only: albedodat, albedice,emisice, emissiv13 USE surfdat_h, only: emisice, emissiv 13 14 USE comgeomfi_h, only: lati 14 15 USE tracer_h, only: noms, rho_co2 … … 32 33 ! ptsrf(ngrid) Surface temperature 33 34 ! 34 ! pdt(ngrid,nlayer) Time derivative before condensation/sublimation of pt35 ! pdt(ngrid,nlayer) Time derivative before condensation/sublimation of pt 35 36 ! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf 36 37 ! pqsurf(ngrid,nq) Sedimentation flux at the surface (kg.m-2.s-1) … … 44 45 ! Both 45 46 ! ---- 46 ! piceco2(ngrid) CO2 ice at the surface (kg/m2)47 ! psolaralb(ngrid) Albedo at the surface48 ! pemisurf(ngrid) Emissivity of the surface47 ! piceco2(ngrid) CO2 ice at the surface (kg/m2) 48 ! albedo(ngrid,L_NSPECTV) Spectral albedo at the surface 49 ! pemisurf(ngrid) Emissivity of the surface 49 50 ! 50 51 ! Authors … … 77 78 REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) 78 79 REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) 80 REAL,INTENT(IN) :: albedo_bareground(ngrid) 81 REAL,INTENT(IN) :: albedo_co2_ice_SPECTV(L_NSPECTV) 79 82 REAL,INTENT(INOUT) :: piceco2(ngrid) 80 REAL,INTENT(OUT) :: psolaralb(ngrid)83 REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV) 81 84 REAL,INTENT(OUT) :: pemisurf(ngrid) 82 85 REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) … … 90 93 ! Local variables 91 94 92 INTEGER l,ig,icap,ilay,it,iq 95 INTEGER l,ig,icap,ilay,it,iq,nw 93 96 94 97 REAL reffrad(ngrid,nlayer) ! radius (m) of the co2 ice particles … … 461 464 piceco2(ig)=0. 462 465 endif 463 if (piceco2(ig).gt.0) then 464 psolaralb(ig) = albedice(icap) 466 if (piceco2(ig) .gt. 1.) then ! CO2 Albedo condition changed to 1 mm coverage. Change by MT2015. 467 DO nw=1,L_NSPECTV 468 albedo(ig,nw) = albedo_co2_ice_SPECTV(nw) 469 ENDDO 465 470 emisref(ig) = emisice(icap) 466 471 else 467 psolaralb(ig) = albedodat(ig) 472 DO nw=1,L_NSPECTV 473 albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine. 474 ENDDO 468 475 emisref(ig) = emissiv 469 476 pemisurf(ig) = emissiv -
trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
r1397 r1482 1 subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf, 1 subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf, & 2 2 qsurf,dqsurf,dqs_hyd,pcapcal, & 3 albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice, & 3 albedo,albedo_bareground, & 4 albedo_snow_SPECTV,albedo_co2_ice_SPECTV, & 5 mu0,pdtsurf,pdtsurf_hyd,hice, & 4 6 pctsrf_sic,sea_ice) 5 7 … … 12 14 USE tracer_h 13 15 use slab_ice_h 14 use callkeys_mod, only: albedosnow,ok_slab_ocean,Tsaldiff,maxicethick,co2cond 16 use callkeys_mod, only: albedosnow,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond 17 use radinc_h, only : L_NSPECTV 15 18 16 19 implicit none … … 25 28 ! ------- 26 29 ! Adapted from LMDTERRE by B. Charnay (2010). Further 27 ! modifications by R. Wordsworth (2010). 30 ! Modifications by R. Wordsworth (2010). 31 ! Spectral albedo by M. Turbet (2015). 28 32 ! 29 33 ! Called by … … 45 49 ! Inputs 46 50 ! ------ 47 real albedoice48 save albedoice49 !$OMP THREADPRIVATE(albedoice)50 51 51 real snowlayer 52 52 parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm … … 72 72 real dqsurf(ngrid,nq), pdtsurf(ngrid) 73 73 real hice(ngrid) 74 real albedo0(ngrid), albedo(ngrid) 74 real albedo(ngrid,L_NSPECTV) 75 real albedo_bareground(ngrid) 76 real albedo_snow_SPECTV(L_NSPECTV) 77 real albedo_co2_ice_SPECTV(L_NSPECTV) 75 78 real pctsrf_sic(ngrid), sea_ice(ngrid) 76 79 … … 85 88 ! ----- 86 89 real a,b,E 87 integer ig,iq, icap ! wld like to remove icap90 integer ig,iq, nw 88 91 real fsnoi, subli, fauxo 89 92 real twater(ngrid) … … 130 133 ! iliq is used in place of igcm_h2o_vap ONLY on the 131 134 ! surface. 132 133 135 ! Soon to be extended to the entire water cycle... 134 135 ! Ice albedo = snow albedo for now136 albedoice=albedosnow137 136 138 137 ! Total ocean surface area … … 188 187 ! albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04) 189 188 ! else 190 albedo(ig) = alb_ocean ! modif Benjamin 189 do nw=1,L_NSPECTV 190 albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in slab_ice_h.F90. Later we could introduce spectral dependency for alb_ocean. 191 enddo 191 192 ! end if 192 193 193 194 194 if(ok_slab_ocean) then195 196 ! Fraction neige (hauteur critique45kg/m2~15cm)197 zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0))198 ! Albedo glace 199 alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) &200 *exp(-sea_ice(ig)/h_alb_ice)201 ! Albedo glace+neige 202 albedo_sic= albedosnow*zfra + alb_ice*(1.0-zfra)203 204 ! Albedo final 205 albedo(ig) = pctsrf_sic(ig)*albedo_sic + (1.-pctsrf_sic(ig))*alb_ocean 206 ! oceanic ice height, just for diagnostics207 hice(ig) = MIN(10.,sea_ice(ig)/rhowater)208 else !ok_slab_ocean195 if(ok_slab_ocean) then 196 197 zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0)) ! Snow Fraction (Critical height 45kg/m2~15cm) 198 alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) & ! Ice Albedo 199 *exp(-sea_ice(ig)/h_alb_ice) 200 ! Albedo final calculation : 201 do nw=1,L_NSPECTV 202 albedo(ig,nw) = pctsrf_sic(ig)* & 203 (albedo_snow_SPECTV(nw)*zfra + alb_ice*(1.0-zfra)) & 204 + (1.-pctsrf_sic(ig))*alb_ocean 205 enddo 206 207 ! Oceanic ice height, just for diagnostics 208 hice(ig) = MIN(10.,sea_ice(ig)/rhowater) 209 else !ok_slab_ocean 209 210 210 211 211 212 ! calculate oceanic ice height including the latent heat of ice formation 212 213 ! hice is the height of oceanic ice with a maximum of maxicethick. 213 hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall 214 215 ! twater(ig) = tsurf(ig) + ptimestep*zdtsurf(ig) & 216 twater(ig) = ztsurf(ig) - hice(ig)*RLFTT*rhowater/pcapcal(ig) 217 ! this is temperature water would have if we melted entire ocean ice layer 218 hicebis(ig) = hice(ig) 219 hice(ig) = 0. 220 221 if(twater(ig) .lt. T_h2O_ice_liq)then 222 E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick) 223 hice(ig) = E/(RLFTT*rhowater) 224 hice(ig) = max(hice(ig),0.0) 225 hice(ig) = min(hice(ig),maxicethick) 226 pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep 227 albedo(ig) = albedoice 228 229 ! if (zqsurf(ig,iice).ge.snowlayer) then 230 ! albedo(ig) = albedoice 231 ! else 232 ! albedo(ig) = albedoocean & 233 ! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer 234 ! endif 235 236 else 237 238 pdtsurf_hyd(ig) = -hicebis(ig)*RLFTT*rhowater/pcapcal(ig)/ptimestep 239 albedo(ig) = alb_ocean 240 241 endif 242 243 zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice)) 244 zqsurf(ig,iice) = hice(ig)*rhowater 245 246 endif!(ok_slab_ocean) 214 hice(ig) = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall 215 216 ! twater(ig) = tsurf(ig) + ptimestep*zdtsurf(ig) & 217 twater(ig) = ztsurf(ig) - hice(ig)*RLFTT*rhowater/pcapcal(ig) 218 ! this is temperature water would have if we melted entire ocean ice layer 219 hicebis(ig) = hice(ig) 220 hice(ig) = 0. 221 222 if(twater(ig) .lt. T_h2O_ice_liq)then 223 E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick) 224 hice(ig) = E/(RLFTT*rhowater) 225 hice(ig) = max(hice(ig),0.0) 226 hice(ig) = min(hice(ig),maxicethick) 227 pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep 228 do nw=1,L_NSPECTV 229 albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015. 230 enddo 231 232 ! if (zqsurf(ig,iice).ge.snowlayer) then 233 ! albedo(ig) = albedoice 234 ! else 235 ! albedo(ig) = albedoocean & 236 ! + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer 237 ! endif 238 239 else 240 241 pdtsurf_hyd(ig) = -hicebis(ig)*RLFTT*rhowater/pcapcal(ig)/ptimestep 242 DO nw=1,L_NSPECTV 243 albedo(ig,nw) = alb_ocean 244 ENDDO 245 246 endif 247 248 zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice)) 249 zqsurf(ig,iice) = hice(ig)*rhowater 250 251 endif!(ok_slab_ocean) 247 252 248 253 … … 299 304 300 305 ! re-calculate continental albedo 301 albedo(ig) = albedo0(ig) ! albedo0 = base values 306 DO nw=1,L_NSPECTV 307 albedo(ig,nw) = albedo_bareground(ig) 308 ENDDO 302 309 if (zqsurf(ig,iice).ge.snowlayer) then 303 albedo(ig) = albedosnow 310 DO nw=1,L_NSPECTV 311 albedo(ig,nw) = albedo_snow_SPECTV(nw) 312 ENDDO 304 313 else 305 albedo(ig) = albedo0(ig) & 306 + (albedosnow - albedo0(ig))*zqsurf(ig,iice)/snowlayer 314 DO nw=1,L_NSPECTV 315 albedo(ig,nw) = albedo_bareground(ig) & 316 + (albedo_snow_SPECTV(nw) - albedo_bareground(ig)) & 317 *zqsurf(ig,iice)/snowlayer 318 ENDDO 307 319 endif 308 320 … … 370 382 if(co2cond)then 371 383 372 icap=1373 384 do ig=1,ngrid 374 if (qsurf(ig,igcm_co2_ice).gt.0) then 375 albedo(ig) = albedice(icap) 376 endif 377 enddo 385 if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015 386 DO nw=1,L_NSPECTV 387 albedo(ig,nw) = albedo_co2_ice_SPECTV(nw) 388 ENDDO 389 endif 390 enddo ! ngrid 378 391 379 endif 392 endif ! co2cond 380 393 381 394 -
trunk/LMDZ.GENERIC/libf/phystd/inifis.F
r1423 r1482 619 619 call getin_p("albedosnow",albedosnow) 620 620 write(*,*) " albedosnow = ",albedosnow 621 622 write(*,*) "CO2 ice albedo ?" 623 albedoco2ice=0.5 ! default value 624 call getin_p("albedoco2ice",albedoco2ice) 625 write(*,*) " albedoco2ice = ",albedoco2ice 621 626 622 627 write(*,*) "Maximum ice thickness ?" -
trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90
r1397 r1482 12 12 use comgeomfi_h, only: area 13 13 use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, & 14 albedice, emisice, emissiv,&14 emisice, emissiv, & 15 15 iceradius, dtemisice, phisfi 16 16 use iostart, only : open_restartphy, close_restartphy, & … … 83 83 84 84 ! Optical properties of polar caps and ground emissivity 85 tab_cntrl(22) = albedice(1) ! Albedo of northern cap ~0.586 tab_cntrl(23) = albedice(2) ! Albedo of southern cap ~0.587 85 tab_cntrl(24) = emisice(1) ! Emissivity of northern cap ~0.95 88 86 tab_cntrl(25) = emisice(2) ! Emissivity of southern cap ~0.95 -
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
r1477 r1482 11 11 use watercommon_h, only : RLVTT, Psat_water,epsi 12 12 use gases_h, only: gnom, gfrac 13 use radcommon_h, only: sigma, glat, grav 13 use radcommon_h, only: sigma, glat, grav, BWNV 14 14 use radii_mod, only: h2o_reffrad, co2_reffrad 15 15 use aerosol_mod, only: iaero_co2, iaero_h2o 16 use surfdat_h, only: phisfi, albedodat,zmea, zstd, zsig, zgam, zthe, &16 use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, & 17 17 dryness, watercaptag 18 18 use comdiurn_h, only: coslat, sinlat, coslon, sinlon … … 190 190 ! ---------------------- 191 191 192 integer day_ini ! Initial date of the run (sol since Ls=0). 193 integer icount ! Counter of calls to physiq during the run. 194 real, dimension(:),allocatable,save :: tsurf ! Surface temperature (K). 195 real, dimension(:,:),allocatable,save :: tsoil ! Sub-surface temperatures (K). 196 real, dimension(:),allocatable,save :: albedo ! Surface albedo. 197 198 !$OMP THREADPRIVATE(tsurf,tsoil,albedo) 199 200 real,dimension(:),allocatable,save :: albedo0 ! Initial surface albedo. 201 real,dimension(:),allocatable,save :: rnat ! Defines the type of the grid (ocean,continent,...). By BC. 202 203 !$OMP THREADPRIVATE(albedo0,rnat) 192 integer day_ini ! Initial date of the run (sol since Ls=0). 193 integer icount ! Counter of calls to physiq during the run. 194 real, dimension(:),allocatable,save :: tsurf ! Surface temperature (K). 195 real, dimension(:,:),allocatable,save :: tsoil ! Sub-surface temperatures (K). 196 real, dimension(:,:),allocatable,save :: albedo ! Surface Spectral albedo. By MT2015. 197 real, dimension(:),allocatable,save :: albedo_equivalent ! Spectral Mean albedo. 198 real, dimension(:),allocatable,save :: albedo_snow_SPECTV ! Snow Spectral albedo. 199 real, dimension(:),allocatable,save :: albedo_co2_ice_SPECTV ! CO2 ice Spectral albedo. 200 201 !$OMP THREADPRIVATE(tsurf,tsoil,albedo,albedo_equivalent,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 202 203 real,dimension(:),allocatable,save :: albedo_bareground ! Bare Ground Albedo. By MT 2015. 204 real,dimension(:),allocatable,save :: rnat ! Defines the type of the grid (ocean,continent,...). By BC. 205 206 !$OMP THREADPRIVATE(albedo_bareground,rnat) 204 207 205 208 real,dimension(:),allocatable,save :: emis ! Thermal IR surface emissivity. … … 228 231 real pw(ngrid,nlayer) ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!) 229 232 230 integer l,ig,ierr,iq 233 integer l,ig,ierr,iq,nw 231 234 232 235 ! FOR DIAGNOSTIC : 233 236 234 real,dimension(:),allocatable,save :: fluxsurf_lw ! Incident Long Wave (IR) surface flux (W.m-2). 235 real,dimension(:),allocatable,save :: fluxsurf_sw ! incident Short Wave (stellar) surface flux (W.m-2). 236 real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2). 237 real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2). 238 real,dimension(:),allocatable,save :: fluxtop_dn ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2). 239 real,dimension(:),allocatable,save :: fluxdyn ! Horizontal heat transport by dynamics (W.m-2). 240 real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)). 241 real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)). 242 real,dimension(:,:),allocatable,save :: zdtlw ! LW heating tendencies (K/s). 243 real,dimension(:,:),allocatable,save :: zdtsw ! SW heating tendencies (K/s). 244 real,dimension(:),allocatable,save :: sensibFlux ! Turbulent flux given by the atmosphere to the surface (W.m-2). 245 246 !$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,& 237 real,dimension(:),allocatable,save :: fluxsurf_lw ! Incident Long Wave (IR) surface flux (W.m-2). 238 real,dimension(:),allocatable,save :: fluxsurf_sw ! Incident Short Wave (stellar) surface flux (W.m-2). 239 real,dimension(:),allocatable,save :: fluxsurfabs_sw ! Absorbed Short Wave (stellar) flux by the surface (W.m-2). 240 real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2). 241 real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2). 242 real,dimension(:),allocatable,save :: fluxtop_dn ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2). 243 real,dimension(:),allocatable,save :: fluxdyn ! Horizontal heat transport by dynamics (W.m-2). 244 real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)). 245 real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)). 246 real,dimension(:,:),allocatable,save :: zdtlw ! LW heating tendencies (K/s). 247 real,dimension(:,:),allocatable,save :: zdtsw ! SW heating tendencies (K/s). 248 real,dimension(:),allocatable,save :: sensibFlux ! Turbulent flux given by the atmosphere to the surface (W.m-2). 249 250 !$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,& 247 251 248 252 !$OMP zdtlw,zdtsw,sensibFlux) … … 376 380 logical clearsky ! For double radiative transfer call. By BC 377 381 378 real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux diagnostics. 379 real tau_col1(ngrid) ! For aerosol optical depth diagnostic. 380 real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV) ! For Outgoing Radiation diagnostics. 382 ! For Clear Sky Case. 383 real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid) ! For SW/LW flux. 384 real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux. 385 real albedo_equivalent1(ngrid) ! For Equivalent albedo calculation. 386 real tau_col1(ngrid) ! For aerosol optical depth diagnostic. 387 real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV) ! For Outgoing Radiation diagnostics. 381 388 real tf, ntf 382 389 … … 435 442 ALLOCATE(tsurf(ngrid)) 436 443 ALLOCATE(tsoil(ngrid,nsoilmx)) 437 ALLOCATE(albedo(ngrid)) 438 ALLOCATE(albedo0(ngrid)) 444 ALLOCATE(albedo(ngrid,L_NSPECTV)) 445 ALLOCATE(albedo_equivalent(ngrid)) 446 ALLOCATE(albedo_snow_SPECTV(L_NSPECTV)) 447 ALLOCATE(albedo_co2_ice_SPECTV(L_NSPECTV)) 448 ALLOCATE(albedo_bareground(ngrid)) 439 449 ALLOCATE(rnat(ngrid)) 440 450 ALLOCATE(emis(ngrid)) … … 457 467 ALLOCATE(fluxsurf_lw(ngrid)) 458 468 ALLOCATE(fluxsurf_sw(ngrid)) 469 ALLOCATE(fluxsurfabs_sw(ngrid)) 459 470 ALLOCATE(fluxtop_lw(ngrid)) 460 471 ALLOCATE(fluxabs_sw(ngrid)) … … 521 532 write (*,*) 'In physiq day_ini =', day_ini 522 533 523 ! Initialize albedo and orbital calculation. 524 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 525 call surfini(ngrid,nq,qsurf,albedo0) 534 ! Initialize albedo calculation. 535 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 536 albedo(:,:)=0.0 537 albedo_bareground(:)=0.0 538 albedo_snow_SPECTV(:)=0.0 539 albedo_co2_ice_SPECTV(:)=0.0 540 call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 541 542 543 ! Initialize orbital calculation. 544 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 526 545 call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) 527 546 528 albedo(:)=albedo0(:)529 547 530 548 if(tlocked)then … … 677 695 if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d. 678 696 call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, & 679 ptimestep,pday+nday,time_phys,area, &680 albedo dat,inertiedat,zmea,zstd,zsig,zgam,zthe)697 ptimestep,pday+nday,time_phys,area, & 698 albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) 681 699 endif 682 700 701 683 702 endif ! end of 'firstcall' 684 703 … … 848 867 ! standard callcorrk 849 868 clearsky=.false. 850 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 851 albedo,emis,mu0,pplev,pplay,pt, & 852 tsurf,fract,dist_star,aerosol,muvar, & 853 zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, & 854 fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, & 855 tau_col,cloudfrac,totcloudfrac, & 856 clearsky,firstcall,lastcall) 857 858 ! Option to call scheme once more for clear regions 869 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 870 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 871 tsurf,fract,dist_star,aerosol,muvar, & 872 zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, & 873 fluxsurfabs_sw,fluxtop_lw, & 874 fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, & 875 tau_col,cloudfrac,totcloudfrac, & 876 clearsky,firstcall,lastcall) 877 878 ! Option to call scheme once more for clear regions 859 879 if(CLFvarying)then 860 880 861 881 ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ... 862 882 clearsky=.true. 863 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 864 albedo,emis,mu0,pplev,pplay,pt, & 865 tsurf,fract,dist_star,aerosol,muvar, & 866 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, & 867 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, & 868 tau_col1,cloudfrac,totcloudfrac, & 883 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 884 albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt, & 885 tsurf,fract,dist_star,aerosol,muvar, & 886 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1, & 887 fluxsurfabs_sw1,fluxtop_lw1, & 888 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1, & 889 tau_col1,cloudfrac,totcloudfrac, & 869 890 clearsky,firstcall,lastcall) 870 891 clearsky = .false. ! just in case. 871 872 892 873 893 ! Sum the fluxes and heating rates from cloudy/clear cases 874 894 do ig=1,ngrid 875 895 tf=totcloudfrac(ig) 876 ntf=1.-tf 877 878 fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig) 879 fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig) 880 fluxtop_lw(ig) = ntf*fluxtop_lw1(ig) + tf*fluxtop_lw(ig) 881 fluxabs_sw(ig) = ntf*fluxabs_sw1(ig) + tf*fluxabs_sw(ig) 882 tau_col(ig) = ntf*tau_col1(ig) + tf*tau_col(ig) 896 ntf=1.-tf 897 fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig) 898 fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig) 899 albedo_equivalent(ig) = ntf*albedo_equivalent1(ig) + tf*albedo_equivalent(ig) 900 fluxsurfabs_sw(ig) = ntf*fluxsurfabs_sw1(ig) + tf*fluxsurfabs_sw(ig) 901 fluxtop_lw(ig) = ntf*fluxtop_lw1(ig) + tf*fluxtop_lw(ig) 902 fluxabs_sw(ig) = ntf*fluxabs_sw1(ig) + tf*fluxabs_sw(ig) 903 tau_col(ig) = ntf*tau_col1(ig) + tf*tau_col(ig) 883 904 884 905 zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer) … … 887 908 OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV) 888 909 OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI) 889 enddo910 enddo 890 911 891 912 endif ! end of CLFvarying. … … 896 917 897 918 898 ! Radiative flux from the sky absorbed by the surface (W.m-2) 919 ! Radiative flux from the sky absorbed by the surface (W.m-2). 899 920 GSR=0.0 900 fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurf _sw(1:ngrid)*(1.-albedo(1:ngrid))921 fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid) 901 922 902 923 !if(noradsurf)then ! no lower surface; SW flux just disappears … … 910 931 911 932 elseif(newtonian)then 912 913 ! II.b Call Newtonian cooling scheme 914 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 933 934 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 935 ! II.b Call Newtonian cooling scheme 936 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 915 937 call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall) 916 938 … … 928 950 endif 929 951 fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid) 930 fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid)) 952 print*,'------------WARNING---WARNING------------' ! by MT2015. 953 print*,'You are in corrk=false mode, ' 954 print*,'and the surface albedo is taken equal to its value at 0.5 micrometers ...' 955 if ( (10000./BWNV(1) .le. 0.5) .or. (10000./BWNV(L_NSPECTV) .ge. 0.5) ) then 956 print*,'0.5 microns band doesnt match your visible bands ! Abort !' 957 call abort 958 endif 959 DO nw=1, L_NSPECTV-1 960 if ( (10000./BWNV(nw) .ge. 0.5) .and. (10000./BWNV(nw+1) .lt. 0.5) ) then 961 fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,nw)) 962 fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid) 963 endif 964 ENDDO 931 965 fluxtop_lw(1:ngrid) = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4 932 966 … … 950 984 call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW) 951 985 call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW) 952 call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk 986 !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk 987 call planetwide_sumval(fluxsurfabs_sw(:)*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk 953 988 call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:)/totarea_planet,dEtotsLW) 954 989 dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:) … … 1166 1201 pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, & 1167 1202 qsurf(1,igcm_co2_ice),albedo,emis, & 1203 albedo_bareground,albedo_co2_ice_SPECTV, & 1168 1204 zdtc,zdtsurfc,pdpsrf,zduc,zdvc, & 1169 1205 zdqc) … … 1461 1497 if(hydrology)then 1462 1498 1463 call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, & 1464 capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, & 1499 call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, & 1500 capcal,albedo,albedo_bareground, & 1501 albedo_snow_SPECTV,albedo_co2_ice_SPECTV, & 1502 mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, & 1465 1503 sea_ice) 1466 1504 … … 1848 1886 call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) 1849 1887 call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) 1850 call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo) 1888 !call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) 1889 !call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) 1851 1890 call wstats(ngrid,"p","Pressure","Pa",3,pplay) 1852 1891 call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) … … 1955 1994 if(callrad.and.(.not.newtonian))then 1956 1995 1957 call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo) 1996 !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) 1997 !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) 1958 1998 call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn) 1959 1999 call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) -
trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F
r961 r1482 24 24 25 25 integer L, NG, NW, NG1,k 26 real*8 rsfv, ubar0, f0pi, btop, bsurf, taumax, eterm 26 real*8 ubar0, f0pi, btop, bsurf, taumax, eterm 27 real*8 rsfv(L_NSPECTV) ! Spectral dependency added by MT2015. 27 28 real*8 FZEROV(L_NSPECTV) 28 29 … … 83 84 84 85 ETERM = MIN(TAUV(L_NLEVRAD,NW,NG),TAUMAX) 85 BSURF = RSFV *UBAR0*STEL(NW)*EXP(-ETERM/UBAR0)86 BSURF = RSFV(NW)*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0) 86 87 87 88 C WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM … … 94 95 95 96 CALL GFLUXV(DTAUV(1,NW,NG),TAUV(1,NW,NG),TAUCUMV(1,NW,NG), 96 * WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV ,97 * WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW), 97 98 * BTOP,BSURF,FMUPV,FMDV,DIFFV,FLUXUP,FLUXDN) 98 99 … … 141 142 142 143 ETERM = MIN(TAUV(L_NLEVRAD,NW,NG),TAUMAX) 143 BSURF = RSFV *UBAR0*STEL(NW)*EXP(-ETERM/UBAR0)144 BSURF = RSFV(NW)*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0) 144 145 145 146 … … 152 153 153 154 CALL GFLUXV(DTAUV(1,NW,NG),TAUV(1,NW,NG),TAUCUMV(1,NW,NG), 154 * WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV ,155 * WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW), 155 156 * BTOP,BSURF,FMUPV,FMDV,DIFFV,FLUXUP,FLUXDN) 156 157 -
trunk/LMDZ.GENERIC/libf/phystd/surfdat_h.F90
r1315 r1482 4 4 implicit none 5 5 6 real,allocatable,dimension(:) :: albedodat ! albedo of bare ground 6 real,allocatable,dimension(:) :: albedodat ! albedo of bare ground stocked in startfi.nc file. 7 7 !$OMP THREADPRIVATE(albedodat) 8 8 ! Ehouarn: moved inertiedat to comsoil.h … … 10 10 real,allocatable,dimension(:) :: phisfi ! geopotential at ground level 11 11 !$OMP THREADPRIVATE(phisfi) 12 real,dimension(2) :: albedice13 12 real,dimension(2) :: emisice ! ice emissivity; 1:Northern hemisphere 2:Southern hemisphere 14 13 real emissiv 15 14 real,dimension(2) :: iceradius, dtemisice 16 !$OMP THREADPRIVATE( albedice,emisice,emissiv,iceradius,dtemisice)15 !$OMP THREADPRIVATE(emisice,emissiv,iceradius,dtemisice) 17 16 real,allocatable,dimension(:) :: zmea,zstd,zsig,zgam,zthe 18 17 !$OMP THREADPRIVATE(zmea,zstd,zsig,zgam,zthe) -
trunk/LMDZ.GENERIC/libf/phystd/surfini.F
r1397 r1482 1 SUBROUTINE surfini(ngrid,nq,qsurf,psolaralb) 1 SUBROUTINE surfini(ngrid,nq,qsurf,albedo,albedo_bareground, 2 & albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 2 3 3 USE surfdat_h, only: albedodat , albedice4 USE surfdat_h, only: albedodat 4 5 USE tracer_h, only: igcm_co2_ice 5 use comgeomfi_h, only: lati6 6 use planetwide_mod, only: planetwide_maxval, planetwide_minval 7 use radinc_h, only : L_NSPECTV 8 use callkeys_mod, only : albedosnow, albedoco2ice 7 9 8 10 IMPLICIT NONE 9 c======================================================================= 10 c 11 c creation des calottes pour l'etat initial 12 c 13 c======================================================================= 14 c----------------------------------------------------------------------- 11 12 13 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 14 cccccccccccccc cccccccccccccc 15 cccccccccccccc Spectral Albedo Initialisation - Routine modified by MT2015. cccccccccccccc 16 cccccccccccccc cccccccccccccc 17 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 18 19 20 c-------------------- 15 21 c Declarations: 16 c ------------- 17 !#include "dimensions.h" 18 !#include "dimphys.h" 19 c 22 c-------------------- 23 20 24 INTEGER,INTENT(IN) :: ngrid 21 25 INTEGER,INTENT(IN) :: nq 22 REAL,INTENT(OUT) :: psolaralb(ngrid) 23 REAL,INTENT(IN) :: qsurf(ngrid,nq) !tracer on surface (kg/m2) 26 REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV) 27 REAL,INTENT(OUT) :: albedo_bareground(ngrid) 28 REAL,INTENT(OUT) :: albedo_snow_SPECTV(L_NSPECTV) 29 REAL,INTENT(OUT) :: albedo_co2_ice_SPECTV(L_NSPECTV) 30 REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg/m2) 24 31 25 INTEGER :: ig, icap32 INTEGER :: ig,nw 26 33 REAL :: min_albedo,max_albedo 27 c 34 28 35 c======================================================================= 29 36 30 31 DO ig=1,ngrid 32 psolaralb(ig)=albedodat(ig) 37 ! Step 1 : Initialisation of the Spectral Albedos. 38 DO nw=1,L_NSPECTV 39 albedo_snow_SPECTV(nw)=albedosnow 40 albedo_co2_ice_SPECTV(nw)=albedoco2ice 33 41 ENDDO 34 42 35 call planetwide_minval(albedodat,min_albedo)36 call planetwide_maxval(albedodat,max_albedo)37 write(*,*) 'surfini: minimum corrected albedo',min_albedo38 write(*,*) 'surfini: maximum corrected albedo',max_albedo39 43 44 ! Step 2 : We get the bare ground albedo from the start files. 45 DO ig=1,ngrid 46 albedo_bareground(ig)=albedodat(ig) 47 DO nw=1,L_NSPECTV 48 albedo(ig,nw)=albedo_bareground(ig) 49 ENDDO 50 ENDDO 51 call planetwide_minval(albedo_bareground,min_albedo) 52 call planetwide_maxval(albedo_bareground,max_albedo) 53 write(*,*) 'surfini: minimum bare ground albedo',min_albedo 54 write(*,*) 'surfini: maximum bare ground albedo',max_albedo 55 56 57 ! Step 3 : We modify the albedo considering some CO2 at the surface. We dont take into account water ice (this is made in hydrol after the first timestep) ... 40 58 if (igcm_co2_ice.ne.0) then 41 ! Change Albedo if there is CO2 ice on the surface 42 DO ig=1,ngrid 43 IF (qsurf(ig,igcm_co2_ice) .GT. 0.) THEN 44 IF(lati(ig).LT.0.) THEN 45 icap=2 ! Southern hemisphere 46 ELSE 47 icap=1 ! Northern hemisphere 48 ENDIF 49 psolaralb(ig) = albedice(icap) 50 END IF 51 ENDDO ! of DO ig=1,ngrid 59 DO ig=1,ngrid 60 IF (qsurf(ig,igcm_co2_ice) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of CO2 ice deposit. 61 DO nw=1,L_NSPECTV 62 albedo(ig,nw)=albedo_co2_ice_SPECTV(nw) 63 ENDDO 64 END IF 65 ENDDO 52 66 else 53 write(*,*) "surfini: No CO2 ice tracer on surface ..." 54 write(*,*) " and therefore no albedo change." 55 endif 67 write(*,*) "surfini: No CO2 ice tracer on surface ..." 68 write(*,*) " and therefore no albedo change." 69 endif 70 call planetwide_minval(albedo,min_albedo) 71 call planetwide_maxval(albedo,max_albedo) 72 write(*,*) 'surfini: minimum corrected initial albedo',min_albedo 73 write(*,*) 'surfini: maximum corrected initial albedo',max_albedo 56 74 57 call planetwide_minval(psolaralb,min_albedo)58 call planetwide_maxval(psolaralb,max_albedo)59 write(*,*) 'surfini: minimum corrected albedo',min_albedo60 write(*,*) 'surfini: maximum corrected albedo',max_albedo61 75 62 76 END -
trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
r1397 r1482 46 46 use ioipsl_getincom , only: getin 47 47 48 use surfdat_h, only: albedice,emisice, iceradius, dtemisice,48 use surfdat_h, only: emisice, iceradius, dtemisice, 49 49 & emissiv 50 50 use comsoil_h, only: volcapa … … 133 133 emin_turb = tab_cntrl(tab0+21) 134 134 c optical properties of polar caps and ground emissivity 135 albedice(1)= tab_cntrl(tab0+22)136 albedice(2)= tab_cntrl(tab0+23)137 135 emisice(1) = tab_cntrl(tab0+24) 138 136 emisice(2) = tab_cntrl(tab0+25) … … 190 188 write(*,5) '(24) emisice(1)',tab_cntrl(tab0+24),emisice(1) 191 189 write(*,5) '(25) emisice(2)',tab_cntrl(tab0+25),emisice(2) 192 write(*,5) '(22) albedice(1)',tab_cntrl(tab0+22),albedice(1)193 write(*,5) '(23) albedice(2)',tab_cntrl(tab0+23),albedice(2)194 190 write(*,6) '(31) iceradius(1)',tab_cntrl(tab0+31),iceradius(1) 195 191 write(*,6) '(32) iceradius(2)',tab_cntrl(tab0+32),iceradius(2) … … 225 221 write(*,*) '(26) emissiv : ground emissivity' 226 222 write(*,*) '(24 et 25) emisice : CO2 ice max emissivity ' 227 write(*,*) '(22 et 23) albedice : CO2 ice cap albedos'228 223 write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow' 229 224 write(*,*) '(33 et 34) dtemisice : time scale for snow', … … 315 310 write(*,*) ' emisice(2) (new value):',emisice(2) 316 311 317 else if (modif(1:len_trim(modif)) .eq. 'albedice') then318 write(*,*) 'current value albedice(1) North:',albedice(1)319 write(*,*) 'enter new value:'320 108 read(*,*,iostat=ierr) albedice(1)321 if(ierr.ne.0) goto 108322 write(*,*)323 write(*,*) ' albedice(1) (new value):',albedice(1)324 write(*,*)325 326 write(*,*) 'current value albedice(2) South:',albedice(2)327 write(*,*) 'enter new value:'328 109 read(*,*,iostat=ierr) albedice(2)329 if(ierr.ne.0) goto 109330 write(*,*)331 write(*,*) ' albedice(2) (new value):',albedice(2)332 333 312 else if (modif(1:len_trim(modif)) .eq. 'iceradius') then 334 313 write(*,*) 'current value iceradius(1) North:',iceradius(1) … … 518 497 write(*,5) '(24) emisice(1)',tab_cntrl(tab0+24),emisice(1) 519 498 write(*,5) '(25) emisice(2)',tab_cntrl(tab0+25),emisice(2) 520 write(*,5) '(22) albedice(1)',tab_cntrl(tab0+22),albedice(1)521 write(*,5) '(23) albedice(2)',tab_cntrl(tab0+23),albedice(2)522 499 write(*,6) '(31) iceradius(1)',tab_cntrl(tab0+31),iceradius(1) 523 500 write(*,6) '(32) iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
Note: See TracChangeset
for help on using the changeset viewer.
