Changeset 1482 for trunk/LMDZ.GENERIC
- Timestamp:
- Oct 14, 2015, 5:05:47 PM (9 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 12 edited
-
README (modified) (1 diff)
-
libf/phystd/callcorrk.F90 (modified) (10 diffs)
-
libf/phystd/callkeys_mod.F90 (modified) (1 diff)
-
libf/phystd/condense_co2.F90 (modified) (6 diffs)
-
libf/phystd/hydrol.F90 (modified) (10 diffs)
-
libf/phystd/inifis.F (modified) (1 diff)
-
libf/phystd/phyredem.F90 (modified) (2 diffs)
-
libf/phystd/physiq.F90 (modified) (18 diffs)
-
libf/phystd/sfluxv.F (modified) (5 diffs)
-
libf/phystd/surfdat_h.F90 (modified) (2 diffs)
-
libf/phystd/surfini.F (modified) (1 diff)
-
libf/phystd/tabfi.F (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r1477 r1482 1085 1085 == 21/09/2015 == MT 1086 1086 - Cleanup of physiq.F90 and changed condense_cloud.F90 to condense_co2.F90 1087 1088 == 14/10/2015 == MT 1089 - Implementation of the Spectral Albedo. Albedo(ngrid) is now Albedo(ngrid,L_NSPECTV) and has a value for each visible band. 1090 - Albedodat/Albedo0 have been removed. We now use Albedo_bareground. 1091 - CO2 Ice Albedo North/South dichotomy has been removed. CO2 ice Albedo is no longer stocked in start file. You can now 1092 prescribe CO2 ice albedo in callphys.def with "albedoco2ice". Its default value is 0.5. -
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.
