#if defined(L08_1A) ! *****************************COPYRIGHT******************************* ! (c) [University of Edinburgh] [2009]. All rights reserved. ! This routine has been licensed to the Met Office for use and ! distribution under the JULES collaboration agreement, subject ! to the terms and conditions set out therein. ! [Met Office Ref SC237] ! *****************************COPYRIGHT******************************* ! SUBROUTINE RELAYERSNOW----------------------------------------------- ! Description: ! Redivide snowpack after changes in depth, conserving mass and energy ! Subroutine Interface: SUBROUTINE relayersnow( land_pts,tile_pts,tile_index, & rgrain0,rho0,sice0,snowfall,snowmass, & tsnow0,nsnow,ds,rgrain,rgrainl,sice, & rho_snow_grnd,sliq,tsnow,rho_snow, & snowdepth ) USE ancil_info, ONLY : & ! imported scalars with intent(in) nsmax ! Maximum possible number of snow layers USE c_0_dg_c, ONLY : & ! imported scalar parameters tm ! temperature at which fresh water freezes and ice melts (K) USE c_perma, ONLY : & ! imported scalar parameters hcapi & ! Specific heat capacity of ice (J/kg/K) ,hcapw ! Specific heat capacity of water (J/kg/K) USE snow_param, ONLY : & ! imported scalars with intent(in) rho_snow_const ! constant density of lying snow (kg per m**3) USE rad_param, ONLY : & ! imported scalars with intent(in) r0 ! Grain size for fresh snow (microns) USE switches, ONLY : & ! imported scalars with intent(in) l_snow_albedo ! Switch for prognostic snow albedo USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook IMPLICIT NONE ! Scalar arguments with intent(in) INTEGER, INTENT(IN) :: & land_pts & ! Total number of land points ,tile_pts ! Number of tile points ! Array arguments with intent(in) INTEGER, INTENT(IN) :: & tile_index(land_pts) ! Index of tile points REAL, INTENT(IN) :: & rgrain0(land_pts) & ! Fresh snow grain size (microns) ,rho0(land_pts) & ! Density of fresh snow (kg/m3) ,sice0(land_pts) & ! Ice content of fresh snow (Kg/m2) ! ! Where NSNOW=0, SICE0 is the mass ! ! of the snowpack. ,snowfall(land_pts) & ! Snowfall reaching the ground (kg/m2) ,snowmass(land_pts) & ! Snow mass on the ground (kg/m2) ,tsnow0(land_pts) ! Temperature of fresh snow (K) ! Array arguments with intent(inout) INTEGER, INTENT(INOUT) :: & nsnow(land_pts) ! Number of snow layers REAL, INTENT(INOUT) :: & ds(land_pts,nsmax) & ! Snow layer thicknesses (m) ,rgrain(land_pts) & ! Snow surface grain size (microns) ,rgrainl(land_pts,nsmax) & ! Snow grain size (microns) ,rho_snow_grnd(land_pts) & ! Snowpack bulk density (kg/m3) ,sice(land_pts,nsmax) & ! Ice content of snow layers (kg/m2) ,sliq(land_pts,nsmax) & ! Liquid content of snow layers (kg/m2) ,tsnow(land_pts,nsmax) ! Snow layer temperatures (K) ! Array arguments with intent(out) REAL, INTENT(OUT) :: & rho_snow(land_pts,nsmax) & ! Snow layer densities (kg/m3) ,snowdepth(land_pts) ! Snow depth (m) ! Local scalar parameters REAL, PARAMETER :: thin_snow_limit = 1.0e-6 ! Maximum snow thickness (m) that is neglected ! during relayering. All contributions ! (mass, energy etc) from that snow are ! neglected. ! Local scalars INTEGER :: & i & ! Land point index ,iznew & ! layer index ,izz & ! layer index ,k & ! Tile point index ,n & ! Snow layer index ,new & ! layer index ,old ! layer index REAL :: & csnow & ! Areal heat capacity of layer (J/K/m2) ,oldremains & ! remaining depth in an old layer (m) ,wt ! weight given to a layer value ! Local arrays INTEGER :: & nold(land_pts) ! Number of layers before adjustment REAL :: & d0(land_pts,0:nsmax) & ! Layer thicknesses before adjustment (m) ! ! D0(:,0) represents new snow if NSNOW>0, ! ! otherwise it is all snow. ,e(0:nsmax) & ! Internal energy before adjustment (J/m2) ,newremains(nsmax) & ! available (unfilled) depth in new layer (m) ,r(0:nsmax) & ! Grain size before adjustment (kg/m2) ,s(0:nsmax) & ! Ice content before adjustment (kg/m2) ,w(0:nsmax) & ! Liquid content before adjustment (kg/m2) ,u(nsmax) ! Layer energy contents INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 REAL(KIND=jprb) :: zhook_handle !----------------------------------------------------------------------- ! External routines used !----------------------------------------------------------------------- EXTERNAL layersnow !----------------------------------------------------------------------- ! Initialise grain variables if they aren't used. !----------------------------------------------------------------------- IF (lhook) CALL dr_hook('RELAYERSNOW',zhook_in,zhook_handle) IF (.NOT.l_snow_albedo) THEN rgrain(:) = r0 rgrainl(:,:) = r0 END IF !----------------------------------------------------------------------- ! Store previous layer thicknesses. !----------------------------------------------------------------------- d0(:,0) = 0.0 nold(:) = nsnow(:) DO k=1,tile_pts i = tile_index(k) IF ( nsnow(i) > 0 ) THEN IF ( sice0(i) > 0.0 ) d0(i,0) = sice0(i) / rho0(i) ELSE IF ( sice0(i) > 0.0 ) THEN ! In this case, SICE0 is the mass of the pack, so use RHO_SNOW_GRND. d0(i,0) = sice0(i) / rho_snow_grnd(i) END IF d0(i,1:nsnow(i)) = ds(i,1:nsnow(i)) END DO !----------------------------------------------------------------------- ! Calculate snowdepth !----------------------------------------------------------------------- DO k=1,tile_pts i = tile_index(k) snowdepth(i) = d0(i,0) DO n=1,nsnow(i) snowdepth(i) = snowdepth(i) + d0(i,n) END DO END DO !----------------------------------------------------------------------- ! Divide snowpack into new layers !----------------------------------------------------------------------- ! DEPENDS ON: layersnow CALL layersnow( land_pts,tile_pts,tile_index,snowdepth,nsnow,ds ) DO k=1,tile_pts i = tile_index(k) IF ( nsnow(i) > 0 ) THEN !----------------------------------------------------------------------- ! Store previous snow layer energy and mass contents !----------------------------------------------------------------------- csnow = sice0(i) * hcapi e(0) = csnow * ( tsnow0(i) - tm ) r(0) = rgrain0(i) ! If the previous timestep had no layers (but could still have snow), ! estimate R(0) using average of previous grain size and that of ! fresh snow. IF ( nold(i) == 0 ) & r(0) = ( 1.0-snowfall(i)/snowmass(i)) * rgrain(i) + & snowfall(i)/snowmass(i) * rgrain0(i) s(0) = sice0(i) w(0) = 0.0 DO n=1,nold(i) csnow = sice(i,n)*hcapi + sliq(i,n)*hcapw e(n) = csnow * ( tsnow(i,n) - tm ) r(n) = rgrainl(i,n) s(n) = sice(i,n) w(n) = sliq(i,n) END DO !----------------------------------------------------------------------- ! Initialise accumulations for new layer values. !----------------------------------------------------------------------- u(:) = 0. sice(i,:) = 0. sliq(i,:) = 0. rgrainl(i,:) = 0. !----------------------------------------------------------------------- ! Set the state of the new layers. !----------------------------------------------------------------------- ! Initialise with all new layers empty. newremains(1:nsnow(i)) = ds(i,1:nsnow(i)) ! Start by filling top new layer. iznew = 1 ! Loop over the old layers. DO old=0,nold(i) ! All of this old layer remains to be reassigned to new layer(s). oldremains = d0(i,old) ! Point to first new layer with remaining space. izz = iznew ! Loop over new layers with remaining space. DO new=izz,nsnow(i) IF ( oldremains > newremains(new) ) THEN !----------------------------------------------------------------------- ! The remaining depth in the new layer will be exhausted by some or ! all of the remaining depth from the old layer. !----------------------------------------------------------------------- ! Decrement old layer by the remaining space in new layer. oldremains = oldremains - newremains(new) ! Add properties from old layer to accumulation for new layer. ! Note that wt is <= 1 since here we have oldRemains>newRemains, ! and oldRemains <= d0. IF ( d0(i,old) > thin_snow_limit ) THEN wt = newremains(new) / d0(i,old) u(new) = u(new) + e(old) * wt sice(i,new) = sice(i,new) + s(old) * wt sliq(i,new) = sliq(i,new) + w(old) * wt rgrainl(i,new) = rgrainl(i,new) + & r(old) * newremains(new) END IF ! Update the pointer to the next new layer with space. izz = new + 1 ELSE !----------------------------------------------------------------------- ! The old layer will be exhausted by this increment. !----------------------------------------------------------------------- ! Decrement available space in the new layer. newremains(new) = newremains(new) - oldremains ! Add properties from old layer to accumulation for new layer. IF ( d0(i,old) > thin_snow_limit ) THEN wt = oldremains / d0(i,old) u(new) = u(new) + e(old) * wt sice(i,new) = sice(i,new) + s(old) * wt sliq(i,new) = sliq(i,new) + w(old) * wt rgrainl(i,new) = rgrainl(i,new) + r(old) * oldremains END IF ! Proceed to the next old layer by exiting from the new layer loop. EXIT END IF END DO ! new layers ! Update pointer to the next new layer with space. iznew = izz END DO ! old layers !----------------------------------------------------------------------- ! Diagnose layer temperatures and densities. !----------------------------------------------------------------------- DO n=1,nsnow(i) csnow = sice(i,n)*hcapi + sliq(i,n)*hcapw tsnow(i,n) = tm + u(n) / csnow rho_snow(i,n) = ( sice(i,n) + sliq(i,n) ) / ds(i,n) rgrainl(i,n) = rgrainl(i,n) / ds(i,n) END DO !----------------------------------------------------------------------- ! Snow surface grain size for radiative calculations !----------------------------------------------------------------------- rgrain(i) = rgrainl(i,1) !----------------------------------------------------------------------- ! Diagnose bulk density of pack. !----------------------------------------------------------------------- rho_snow_grnd(i) = snowmass(i) / snowdepth(i) ELSE !----------------------------------------------------------------------- ! Set bulk density of pack to a constant value if there is (effectively) ! no snow. This density is then used the next time a shallow pack forms. ! This is necessary if nsmax>0, but redundant if nsmax=0 (constant ! density). We could also recalculate snowdepth, to make it consistent ! with the revised density, but we're not bothering as depths are tiny! !----------------------------------------------------------------------- IF ( snowmass(i) < 1.0e-9 ) rho_snow_grnd(i) = rho_snow_const END IF ! nsnow !----------------------------------------------------------------------- ! Set values for unused snow layers. ! Note: not needed for algorithm, but clearer to follow. !----------------------------------------------------------------------- IF ( nsnow(i) < nsmax ) THEN n = nsnow(i) + 1 rgrainl(i,n:) = r0 rho_snow(i,n:) = 0.0 sice(i,n:) = 0. sliq(i,n:) = 0. tsnow(i,n:) = tm END IF END DO ! K (points) IF (lhook) CALL dr_hook('RELAYERSNOW',zhook_out,zhook_handle) RETURN END SUBROUTINE relayersnow #endif