Changeset 1482 for trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
- Timestamp:
- Oct 14, 2015, 5:05:47 PM (9 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.
