- Timestamp:
- Nov 19, 2021, 4:58:59 PM (3 years ago)
- Location:
- LMDZ6/branches/Ocean_skin
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/inlandsis/surf_inlandsis_mod.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Ocean_skin
- Property svn:mergeinfo changed
-
LMDZ6/branches/Ocean_skin/libf/phylmd/inlandsis/surf_inlandsis_mod.F90
r3792 r4013 1 1 MODULE surf_inlandsis_mod 2 2 3 IMPLICIT NONE 4 3 IMPLICIT NONE 4 5 CONTAINS 6 7 8 SUBROUTINE surf_inlandsis(knon, rlon, rlat, ikl2i, itime, dtime, debut, lafin, & 9 rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, & 10 precip_rain, precip_snow, & 11 zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, & 12 rugos, snow_cont_air, alb_soil, alt, slope, cloudf, & 13 radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, & 14 AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, & 15 runoff_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s,dflux_l, & 16 tsurf_new, alb1, alb2, alb3, alb6, emis_new, z0m, z0h, qsurf) 17 18 ! | | 19 ! | SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow| 20 ! | (INLANDSIS) | 21 ! | SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme) | 22 ! | surface scheme of the Modele Atmospherique Regional (MAR) | 23 ! | Author: Heinz Juergen Punge, LSCE June 2009 | 24 ! | based on the MAR-SISVAT interface by Hubert Gallee | 25 ! | Updated by Etienne Vignon, Cecile Agosta | 26 ! | | 27 ! +------------------------------------------------------------------------+ 28 ! | 29 ! | In the current setup, SISVAT is used only to model the land ice | 30 ! | part of the surface; hence it is called with the compressed variables| 31 ! | from pbl_surface, and only by the surf_landice routine. | 32 ! | | 33 ! | In this interface it is assumed that the partitioning of the soil, | 34 ! | and hence the number of grid points is constant during a simulation, | 35 ! | hence eg. snow properties remain stored in the global SISVAT | 36 ! | variables between the calls and don't need to be handed over as | 37 ! | arguments. When the partitioning is supposed to change, make sure to | 38 ! | update the variables. | 39 ! | | 40 ! | INPUT (via MODULES VARxSV, VARySV, VARtSV ...) | 41 ! | ^^^^^ xxxxSV: SISVAT/LMDZ interfacing variables | 42 ! | | 43 ! +------------------------------------------------------------------------+ 44 45 USE dimphy 46 USE VAR_SV 47 USE VARdSV 48 USE VARxSV 49 USE VARySV 50 USE VARtSV 51 USE VARphy 52 USE surface_data, only : iflag_tsurf_inlandsis, SnoMod, BloMod, ok_outfor 53 54 IMPLICIT NONE 55 56 ! +--INTERFACE Variables 57 ! + =================== 58 ! include "dimsoil.h" 59 60 ! +--Global Variables 61 ! + ================ 62 ! Input Variables for SISVAT 63 INTEGER, INTENT(IN) :: knon 64 INTEGER, INTENT(IN) :: itime 65 REAL, INTENT(IN) :: dtime 66 LOGICAL, INTENT(IN) :: debut ! true if first step 67 LOGICAL, INTENT(IN) :: lafin ! true if last step 68 69 INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i ! Index Decompression 70 REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat 71 REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cos sol. zenith angle 72 REAL, DIMENSION(klon), INTENT(IN) :: swdown ! 73 REAL, DIMENSION(klon), INTENT(IN) :: lwdown ! 74 REAL, DIMENSION(klon), INTENT(IN) :: albedo_old 75 REAL, DIMENSION(klon), INTENT(IN) :: pexner ! Exner potential 76 REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow 77 REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo 78 REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps, p1lay 79 REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf 80 REAL, DIMENSION(klon), INTENT(IN) :: rugos 81 REAL, DIMENSION(klon), INTENT(IN) :: snow_cont_air 82 REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope 83 REAL, DIMENSION(klon), INTENT(IN) :: alt ! surface elevation 84 REAL, DIMENSION(klon), INTENT(IN) :: cloudf 85 REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ 86 REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ 87 REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh 88 REAL, DIMENSION(klon), INTENT(IN) :: ustar ! friction velocity 89 90 ! Variables exchanged between LMDZ and SISVAT 91 REAL, DIMENSION(klon), INTENT(IN) :: radsol ! Surface absorbed rad. 92 REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! Tot snow mass [kg/m2] 93 REAL, DIMENSION(klon), INTENT(INOUT) :: zfra ! snwo surface fraction [0-1] 94 REAL, DIMENSION(klon, nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature 95 REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! Soil Water Content 96 REAL, DIMENSION(klon), INTENT(INOUT) :: z0m ! Momentum Roughn Lgt 97 REAL, DIMENSION(klon), INTENT(INOUT) :: z0h ! Momentum Roughn Lgt 98 99 ! Output Variables for LMDZ 100 REAL, DIMENSION(klon), INTENT(OUT) :: alb1 ! Albedo SW 101 REAL, DIMENSION(klon), INTENT(OUT) :: alb2, alb3 ! Albedo NIR and LW 102 REAL, DIMENSION(klon,6), INTENT(OUT) :: alb6 ! 6 band Albedo 103 REAL, DIMENSION(klon), INTENT(OUT) :: emis_new ! Surface Emissivity 104 REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff 105 REAL, DIMENSION(klon), INTENT(OUT) :: ffonte ! enthalpy flux due to surface melting 106 REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte ! water flux due to surface melting 107 REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s ! d/dT sens. ht flux 108 REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l ! d/dT latent ht flux 109 REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens ! Sensible ht flux 110 REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat ! Latent heat flux 111 REAL, DIMENSION(klon), INTENT(OUT) :: evap ! Evaporation 112 REAL, DIMENSION(klon), INTENT(OUT) :: erod ! Erosion of surface snow (flux) 113 REAL, DIMENSION(klon), INTENT(OUT) :: agesno ! Snow age (top layer) 114 REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature 115 REAL, DIMENSION(klon), INTENT(OUT) :: qsurf ! Surface Humidity 116 117 ! Specific INLANDIS outputs 118 REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! Total H2O snow[kg/m2] 119 REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! Snow height (m) 120 REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! Snow passed to ice 121 REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! Snow in model (kg/m2) 122 123 ! +--Internal Variables 124 ! + =================== 125 126 CHARACTER(len = 20) :: fn_outfor ! Name for output file 127 CHARACTER (len = 80) :: abort_message 128 CHARACTER (len = 20) :: modname = 'surf_inlandsis_mod' 129 130 INTEGER :: i, ig, ikl, isl, isn, nt 131 INTEGER :: gp_outfor, un_outfor 132 REAL, PARAMETER :: f1 = 0.5 133 REAL, PARAMETER :: sn_upp = 10000., sn_low = 500. 134 REAL, PARAMETER :: sn_add = 400., sn_div = 2. 135 ! snow mass upper,lower limit, 136 ! added mass/division lowest layer 137 REAL, PARAMETER :: c1_zuo = 12.960e+4, c2_zuo = 2.160e+6 138 REAL, PARAMETER :: c3_zuo = 1.400e+2, czemin = 1.e-3 139 ! Parameters for drainage 140 ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ ! Tuning 141 ! +... Run Off Parameters 142 ! + 86400*1.5 day ...*25 days (Modif. ETH Camp: 86400*0.3day) 143 ! + (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317) 144 145 REAL, DIMENSION(klon) :: eps0SL ! surface Emissivity 146 REAL :: zsigma, Ua_min, Us_min, lati 147 REAL, PARAMETER :: cdmax=0.05 148 REAL :: lambda ! Par. soil discret. 149 REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2 ! Soil layer thicknesses 150 !$OMP THREADPRIVATE(dz1,dz2) 151 LOGICAL, SAVE :: firstcall 152 !$OMP THREADPRIVATE(firstcall) 153 154 INTEGER :: iso 155 LOGICAL :: file_exists 156 CHARACTER(len = 20) :: fichnom 157 LOGICAL :: is_init_domec 158 ! CA initialization 159 ! dz_profil_15 : 1 m in 15 layers [m] 160 real, parameter :: dz_profil_15(15) = (/0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, & 161 0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.17/) 162 ! mean_temp : mean annual surface temperature [K] 163 real, dimension(klon) :: mean_temp 164 ! mean_dens : mean surface density [kg/m3] 165 real, dimension(klon) :: mean_dens 166 ! lat_scale : temperature lapse rate against latitude [K degree-1] 167 real :: lat_scale 168 ! sh_scale : temperature lapse rate against altitude [K km-1] 169 real :: sh_scale 170 ! variables for density profile 171 ! E0, E1 : exponent 172 real :: E0, E1 173 ! depth at which 550 kg m-3 is reached [m] 174 real :: z550 175 ! depths of snow layers 176 real :: depth, snow_depth, distup 177 ! number of initial snow layers 178 integer :: nb_snow_layer 179 ! For density calc. 180 real :: alpha0, alpha1, ln_smb 181 ! theoritical densities [kg m-3] 182 real :: rho0, rho1, rho1_550 183 ! constants for density profile 184 ! C0, C1 : constant, 0.07 for z <= 550 kg m-3 185 real, parameter :: C0 = 0.07 186 real, parameter :: C1 = 0.03 187 ! rho_i : ice density [kg m-3] 188 real, parameter :: rho_ice = 917. 189 ! E_c : activation energy [J mol-1] 190 real, parameter :: E_c = 60000. 191 ! E_g : activation energy [J mol-1] 192 real, parameter :: E_g = 42400. 193 ! R : gas constant [J mol-1 K-1] 194 real, parameter :: R = 8.3144621 195 196 197 198 199 200 ! + PROGRAM START 201 ! + ----------------------------------------- 202 203 zsigma = 1000. 204 dt__SV = dtime 205 206 IF (debut) THEN 207 firstcall = .TRUE. 208 INI_SV = .false. 209 ELSE 210 firstcall = .false. 211 INI_SV = .true. 212 END IF 213 214 IF (ok_outfor) THEN 215 un_outfor = 51 ! unit number for point output file 216 gp_outfor = 1 ! grid point number for point output 1 for 1D, 273 for zoom-nudg DC 217 fn_outfor = 'outfor_SV.dat' 218 END IF 219 220 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 221 ! + INITIALISATION: BEGIN +++ 222 ! + ----------------------------------------- 223 IF (firstcall) THEN 224 225 ! +--Array size 226 ! + ----------------------- 227 228 klonv = klon 229 knonv = knon 230 write(*, *) 'ikl, lon and lat in INLANDSIS' 231 232 DO ikl = 1, knon 233 i=ikl2i(ikl) 234 write(*, *) 'ikl=', ikl, 'rlon=', rlon(i), 'rlat=', rlat(i) 235 END DO 236 237 ! +--Variables initizialisation 238 ! + --------------------------- 239 240 CALL INIT_VARtSV 241 CALL INIT_VARxSV 242 CALL INIT_VARySV 243 244 245 246 ! +--Surface Fall Line Slope 247 ! + ----------------------- 248 IF (SnoMod) THEN 249 DO ikl = 1, knon 250 slopSV(ikl) = slope(ikl) 251 SWf_SV(ikl) = & ! Normalized Decay of the 252 exp(-dt__SV & ! Surficial Water Content 253 / (c1_zuo & !(Zuo and Oerlemans 1996, 254 + c2_zuo * exp(-c3_zuo * abs(slopSV(ikl))))) ! J.Glacio. 42, 305--317) 255 END DO 256 END IF 257 258 259 260 ! +--Soil layer thickness . Compute soil discretization (as for LMDZ) 261 ! + ---------------------------------------------------------------- 262 ! write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx 263 CALL get_soil_levels(dz1, dz2, lambda) 264 265 lambSV = lambda 266 dz1_SV(1:knon, 1:) = 0. 267 dz2_SV(1:knon, 1:) = 0. 268 269 DO isl = -nsol, 0 270 dz_dSV(isl) = 0.5e-3 * dz2(1 - isl) ! Soil layer thickness 271 DO ikl = 1, knon 272 dz1_SV(ikl, isl) = dz1(1 - isl) !1.e-3* 273 dz2_SV(ikl, isl) = dz2(1 - isl) !1.e-3* 274 END DO 275 END DO 276 277 278 ! Set variables 279 ! ============= 280 DO ikl = 1, knon 281 ! LSmask : Land/Sea Mask 282 LSmask(ikl) = 1 283 ! isotSV : Soil Type -> 12 = ice 284 isotSV(ikl) = 12 285 ! iWaFSV : Soil Drainage (1,0)=(y,n) 286 iWaFSV(ikl) = 1 287 ! eps0SL : Surface Emissivity 288 eps0SL(ikl) = 1. 289 ! alb0SV : Soil Albedo 290 alb0SV(ikl) = alb_soil(ikl) 291 ! Tsf_SV : Surface Temperature, must be bellow freezing 292 Tsf_SV(ikl) = min(temp_air(ikl), TfSnow) 293 END DO 294 295 ! +--Initialization of soil and snow variables in case startsis is not read 296 ! + ---------------------------------------------------------------------- 297 298 is_init_domec=.FALSE. 299 300 301 IF (is_init_domec) THEN 302 ! Coarse initilization inspired from vertcical profiles at Dome C, 303 ! Antarctic Plateaui (10m of snow, 19 levels) 304 305 DO ikl = 1,knon 306 ! + Soil 307 DO isl = -nsol,0 308 TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2) !temp_air(ikl) 309 !tsoil(ikl,1-isl) Soil Temperature 310 !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2) 311 eta_SV(ikl,isl) = epsi !etasoil(ikl,1-isl) Soil Water[m3/m3] 312 ro__SV(ikl,isl) = rhoIce !rosoil(ikl,1-isl) volumic mass 313 END DO 314 315 316 ! Snow 317 isnoSV(ikl) = 19 318 istoSV(ikl, 1:isnoSV(ikl)) = 100 319 ro__SV(ikl, 1:isnoSV(ikl)) = 350. 320 eta_SV(ikl, 1:isnoSV(ikl)) = epsi 321 TsisSV(ikl, 1:isnoSV(ikl)) = min(tsoil(ikl, 1), TfSnow - 0.2) 322 G1snSV(ikl, 1:isnoSV(ikl)) = 99. 323 G2snSV(ikl, 1:isnoSV(ikl)) = 2. 324 agsnSV(ikl, 1:isnoSV(ikl)) = 50. 325 dzsnSV(ikl, 19) = 0.015 326 dzsnSV(ikl, 18) = 0.015 327 dzsnSV(ikl, 17) = 0.020 328 dzsnSV(ikl, 16) = 0.030 329 dzsnSV(ikl, 15) = 0.040 330 dzsnSV(ikl, 14) = 0.060 331 dzsnSV(ikl, 13) = 0.080 332 dzsnSV(ikl, 12) = 0.110 333 dzsnSV(ikl, 11) = 0.150 334 dzsnSV(ikl, 10) = 0.200 335 dzsnSV(ikl, 9) = 0.300 336 dzsnSV(ikl, 8) = 0.420 337 dzsnSV(ikl, 7) = 0.780 338 dzsnSV(ikl, 6) = 1.020 339 dzsnSV(ikl, 5) = 0.980 340 dzsnSV(ikl, 4) = 1.020 341 dzsnSV(ikl, 3) = 3.970 342 dzsnSV(ikl, 2) = 1.020 343 dzsnSV(ikl, 1) = 1.020 344 345 END DO 346 ELSE 347 348 ! Initilialisation with climatological temperature and density 349 ! profiles as in MAR. Methodology developed by Cecile Agosta 350 351 ! initialize with 0., for unused snow layers 352 dzsnSV = 0. 353 G1snSV = 0. 354 G2snSV = 0. 355 istoSV = 0 356 TsisSV = 0. 357 358 359 ! initialize mean variables (unrealistic) 360 mean_temp = TfSnow 361 mean_dens = 300. 362 ! loop on grid cells 363 DO ikl = 1, knon 364 lati=rlat(ikl2i(ikl)) 365 ! approximations for mean_temp and mean_dens 366 ! from Feulner et al., 2013 (DOI: 10.1175/JCLI-D-12-00636.1) 367 ! Fig. 3 and 5 : the lapse rate vs. latitude at high latitude is about 0.55 °C °lat-1 368 ! with a moist-adiabatic lapse rate of 5 °C km-1 everywhere except for Antarctica, 369 ! for Antarctica, a dry-adiabatic lapse rate of 9.8 °C km-1 is assumed. 370 if (lati > 60.) then 371 ! CA todo : add longitude bounds 372 ! Greenland mean temperature : function of altitude and latitude 373 ! for altitudes 0. to 1000. m, lat_scale varies from 0.9 to 0.75 °C °lat-1 374 lat_scale = (0.75 - 0.9) / 1000. * alt(ikl) + 0.9 375 lat_scale = max(min(lat_scale, 0.9), 0.75) 376 ! sh_scale equals the environmental lapse rate : 6.5 °C km-1 377 sh_scale = 6.5 378 mean_temp(ikl) = TfSnow + 1.5 - sh_scale * alt(ikl) / 1000. - lat_scale * (lati - 60.) 379 ! surface density: Fausto et al. 2018, https://doi.org/10.3389/feart.2018.00051 380 mean_dens(ikl) = 315. 381 else if (lati < -60.) then 382 ! Antarctica mean temperature : function of altitude and latitude 383 ! for altitudes 0. to 500. m, lat_scale varies from 1.3 to 0.6 °C °lat-1 384 lat_scale = (0.6 - 1.3) / 500. * alt(ikl) + 1.3 385 lat_scale = max(min(lat_scale, 1.3), 0.6) 386 ! for altitudes 0. to 500. m, sh_scale varies from 6.5 to 9.8 °C km-1 387 sh_scale = (9.8 - 6.5) / 500. * alt(ikl) + 6.5 388 sh_scale = max(min(sh_scale, 9.8), 6.5) 389 mean_temp(ikl) = TfSnow - 7. - sh_scale * alt(ikl) / 1000. + lat_scale * (lati + 60.) 390 ! Antarctica surface density : function of mean annual temperature 391 ! surface density of 350. kg m-3 at Dome C and 450. kg m-3 at Prud'homme (Agosta et al. 2013) 392 ! 350 kg m-3 is a typical value for the Antarctic plateau around 3200 m. 393 ! Weinhart et al 2020 https://doi.org/10.5194/tc-14-3663-2020 and Sugiyama et al. 2011 oi: 10.3189/2012JoG11J201 394 ! 320 kg m-3 is reached at Dome A, 4100 m a.s.l. 395 ! Dome C : st_ant_param(3233, -75.1) = -47.7 396 ! Dumont d'Urville : st_ant_param(0, -66.66) = -15.7 397 mean_dens(ikl) = (450. - 320.) / (-15.7 + 47.7) * (mean_temp(ikl) - TfSnow + 15.7) + 450. 398 mean_dens(ikl) = min(450., max(320., mean_dens(ikl))) 399 else 400 401 ! write(*, *) 'Attention: temperature initialization is only defined for Greenland and Antarctica' 402 403 mean_dens(ikl) =350. 404 mean_temp(ikl) = min(tsoil(ikl,1),TfSnow-0.2) 405 406 !abort_message='temperature initialization is only defined for Greenland and Antarctica' 407 !CALL abort_physic(modname,abort_message,1) 408 409 end if 5 410 6 CONTAINS 7 8 9 10 SUBROUTINE surf_inlandsis(knon,rlon,rlat, ikl2i, itime, dtime, debut, lafin, & 11 rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, & 12 precip_rain, precip_snow, precip_snow_adv, snow_adv, & 13 zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, & 14 rugos, snow_cont_air, alb_soil, slope, cloudf, & 15 radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, & 16 AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, & 17 runoff_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, & 18 tsurf_new, alb1, alb2, alb3, & 19 emis_new, z0m, z0h, qsurf) 20 21 ! +------------------------------------------------------------------------+ 22 ! | | 23 ! | SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow| 24 ! | (INLANDSIS) | 25 ! | SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme) | 26 ! | surface scheme of the Modele Atmospherique Regional (MAR) | 27 ! | Author: Heinz Juergen Punge, LSCE June 2009 | 28 ! | based on the MAR-SISVAT interface by Hubert Gallee | 29 ! | Update Etienne Vignon, LMD, Novembre 2020 | 30 ! | | 31 ! +------------------------------------------------------------------------+ 32 ! | 33 ! | In the current setup, SISVAT is used only to model the land ice | 34 ! | part of the surface; hence it is called with the compressed variables| 35 ! | from pbl_surface, and only by the surf_landice routine. | 36 ! | | 37 ! | In this interface it is assumed that the partitioning of the soil, | 38 ! | and hence the number of grid points is constant during a simulation, | 39 ! | hence eg. snow properties remain stored in the global SISVAT | 40 ! | variables between the calls and don't need to be handed over as | 41 ! | arguments. When the partitioning is supposed to change, make sure to | 42 ! | update the variables. | 43 ! | | 44 ! | INPUT | 45 ! | SnoMod: Snow Pack is set up when .T. | 46 ! | reaLBC: Update Bound.Condit.when .T. | 47 ! | | 48 ! | INPUT (via MODULES VARxSV, VARySV, VARtSV) | 49 ! | ^^^^^ xxxxSV: SISVAT/LMDZ interfacing variables | 50 ! | | 51 ! | Preprocessing Option: SISVAT PHYSICS | 52 ! | ^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^ | 53 ! | # #HY | 54 ! | # #SN: Snow Model | 55 ! | # #BS: Blowing Snow Parameterization | 56 ! +------------------------------------------------------------------------+ 57 58 USE dimphy 59 USE VAR_SV 60 USE VARdSV 61 USE VARxSV 62 USE VARySV 63 USE VARtSV 64 USE VARphy 65 USE surface_data, only: iflag_tsurf_inlandsis,SnoMod,BloMod,ok_outfor 66 67 IMPLICIT NONE 68 69 ! +--INTERFACE Variables 70 ! + =================== 71 72 ! include "dimsoil.h" 73 74 75 ! +--Global Variables 76 ! + ================ 77 ! Input Variables for SISVAT 78 INTEGER, INTENT(IN) :: knon 79 INTEGER, INTENT(IN) :: itime 80 REAL, INTENT(IN) :: dtime 81 LOGICAL, INTENT(IN) :: debut ! true if first step 82 LOGICAL, INTENT(IN) :: lafin ! true if last step 83 84 INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i ! Index Decompression 85 REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat 86 REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cos sol. zenith angle 87 REAL, DIMENSION(klon), INTENT(IN) :: swdown ! 88 REAL, DIMENSION(klon), INTENT(IN) :: lwdown ! 89 REAL, DIMENSION(klon), INTENT(IN) :: albedo_old 90 REAL, DIMENSION(klon), INTENT(IN) :: pexner ! Exner potential 91 REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow 92 REAL, DIMENSION(klon), INTENT(IN) :: precip_snow_adv, snow_adv 93 !Snow Drift 94 REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo 95 REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps,p1lay 96 REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf 97 REAL, DIMENSION(klon), INTENT(IN) :: rugos,snow_cont_air 98 REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope 99 REAL, DIMENSION(klon), INTENT(IN) :: cloudf 100 REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ 101 REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ 102 REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh 103 REAL, DIMENSION(klon), INTENT(IN) :: ustar ! friction velocity 104 105 ! Variables exchanged between LMDZ and SISVAT 106 REAL, DIMENSION(klon), INTENT(IN) :: radsol ! Surface absorbed rad. 107 REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! Tot snow mass [kg/m2] 108 REAL, DIMENSION(klon), INTENT(INOUT) :: zfra ! snwo surface fraction [0-1] 109 REAL, DIMENSION(klon,nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature 110 REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! Soil Water Content 111 REAL, DIMENSION(klon), INTENT(INOUT) :: z0m ! Momentum Roughn Lgt 112 REAL, DIMENSION(klon), INTENT(INOUT) :: z0h ! Momentum Roughn Lgt 113 114 115 ! Output Variables for LMDZ 116 REAL, DIMENSION(klon), INTENT(OUT) :: alb1 ! Albedo SW 117 REAL, DIMENSION(klon), INTENT(OUT) :: alb2,alb3 ! Albedo NIR and LW 118 REAL, DIMENSION(klon), INTENT(OUT) :: emis_new ! Surface Emissivity 119 REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff 120 REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s ! d/dT sens. ht flux 121 REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l ! d/dT latent ht flux 122 REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens ! Sensible ht flux 123 REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat ! Latent heat flux 124 REAL, DIMENSION(klon), INTENT(OUT) :: evap ! Evaporation 125 REAL, DIMENSION(klon), INTENT(OUT) :: agesno ! Snow age (top layer) 126 REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature 127 REAL, DIMENSION(klon), INTENT(OUT) :: qsurf ! Surface Humidity 128 129 ! Specific INLANDIS outputs 130 131 REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! Total H2O snow[kg/m2] 132 REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! Snow height (m) 133 REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! Snow passed to ice 134 REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! Snow in model (kg/m2) 135 136 137 138 139 ! +--Internal Variables 140 ! + =================== 141 142 CHARACTER(len=20) :: fn_outfor ! Name for output file 143 INTEGER :: i, ig, ikl, isl, isn, nt 144 INTEGER :: gp_outfor, un_outfor 145 REAL, PARAMETER :: f1=0.5 146 REAL, PARAMETER :: sn_upp=5000.,sn_low=500. 147 REAL, PARAMETER :: sn_add=400.,sn_div=2. 148 ! snow mass upper,lower limit, 149 ! added mass/division lowest layer 150 REAL, PARAMETER :: c1_zuo=12.960e+4, c2_zuo=2.160e+6 151 REAL, PARAMETER :: c3_zuo=1.400e+2, czemin=1.e-3 152 ! Parameters for drainage 153 ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ ! Tuning 154 ! +... Run Off Parameters 155 ! + 86400*1.5 day ...*25 days (Modif. ETH Camp: 86400*0.3day) 156 ! + (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317) 157 158 REAL, DIMENSION(klon) :: eps0SL ! surface Emissivity 159 REAL :: zsigma, Ua_min, Us_min 160 REAL :: lambda ! Par. soil discret. 161 REAL, DIMENSION(nsoilmx), SAVE :: dz1,dz2 ! Soil layer thicknesses 162 !$OMP THREADPRIVATE(dz1,dz2) 163 LOGICAL, SAVE :: firstcall 164 !$OMP THREADPRIVATE(firstcall) 165 166 167 168 ! +--Internal Variables 169 ! + ================== 170 171 INTEGER :: iso 172 LOGICAL :: file_exists 173 CHARACTER(len=20) :: fichnom 174 !======================================================================== 175 176 PRINT*, 'je rentre dans inlandsis' 177 178 zsigma=1000. 179 dt__SV=dtime 180 181 182 183 ! write(*,*)'Start of simulation? ',debut !hj 184 185 IF (debut) THEN 186 firstcall=.TRUE. 187 INI_SV=.false. 188 189 ELSE 190 firstcall=.false. 191 INI_SV=.true. 192 END IF 193 194 195 196 197 IF (ok_outfor) THEN 198 un_outfor=51 ! unit number for point output file 199 gp_outfor= 1 ! grid point number for point output 200 fn_outfor='outfor_SV.dat' 201 END IF 202 203 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 204 205 206 207 208 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 209 ! + INITIALISATION: BEGIN +++ 210 ! + ------------------------- 211 ! + 212 ! + Compute soil discretization (as for LMDZ) 213 ! + ----------------------------------------- 214 IF (firstcall) THEN 215 216 ! +--Array size 217 klonv=klon 218 knonv=knon 219 220 221 write(*,*)'klon',klon,'klonv',klonv,'knon',knon,'nsol',nsol,'nsno',nsno 222 223 224 CALL INIT_VARtSV 225 CALL INIT_VARxSV 226 CALL INIT_VARySV 227 228 eps0SL(:)=0. 229 230 231 ! +--Soil layer thickness 232 ! + ----------------------- 233 ! write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx 234 CALL get_soil_levels(dz1,dz2,lambda) 235 236 237 lambSV=lambda 238 dz1_SV(1:knon,1:) = 0. 239 dz2_SV(1:knon,1:) = 0. 240 241 DO isl = -nsol,0 242 dz_dSV(isl) = 0.5e-3*dz2(1-isl) ! Soil layer thickness 243 DO ikl=1,knon 244 dz1_SV(ikl,isl) = dz1(1-isl) !1.e-3* 245 dz2_SV(ikl,isl) = dz2(1-isl) !1.e-3* 246 END DO 411 ! mean_temp is defined for ice ground only 412 mean_temp(ikl) = min(mean_temp(ikl), TfSnow - 0.2) 413 414 ! Soil layers 415 ! =========== 416 DO isl = -nsol, 0 417 ! TsisSV : Temperature [K] 418 TsisSV(ikl, isl) = mean_temp(ikl) 419 ! eta_SV : Soil Water [m3/m3] 420 eta_SV(ikl, isl) = epsi 421 ! ro__SV : Volumic Mass [kg/m3] 422 ro__SV(ikl, isl) = rhoIce 423 END DO 424 425 ! Snow layers 426 ! =========== 427 ! snow_depth : initial snow depth 428 snow_depth = 20. 429 ! nb_snow_layer : initial nb of snow layers 430 nb_snow_layer = 15 431 ! isnoSV : total nb of snow layers 432 isnoSV(ikl) = nb_snow_layer 433 ! depth : depth of each layer 434 depth = snow_depth 435 do isl = 1, nb_snow_layer 436 ! dzsnSV : snow layer thickness 437 dzsnSV(ikl, isl) = max(0.01, snow_depth * dz_profil_15(nb_snow_layer - isl + 1)) 438 ! G1snSV : dendricity (<0) or sphericity (>0) : 99. = sperical 439 G1snSV(ikl, isl) = 99. 440 ! G2snSV : Sphericity (>0) or Size [1/10 mm] : 2. = small grain size 441 G2snSV(ikl, isl) = 3. 442 agsnSV(ikl, isl) = 0. 443 istoSV(ikl, isl) = 0 444 ! eta_SV : Liquid Water Content [m3/m3] 445 eta_SV(ikl, isl) = 0. 446 ! distance to surface 447 depth = depth - dzsnSV(ikl,isl) / 2. 448 distup = min(1., max(0., depth / snow_depth)) 449 ! TsisSV : Temperature [K], square interpolation between Tsf_SV (surface) and mean_temp (bottom) 450 TsisSV(ikl, isl) = Tsf_SV(ikl) * (1. - distup**2) + mean_temp(ikl) * distup**2 451 ! firn density : densification formulas from : 452 ! Ligtenberg et al 2011 eq. (6) (www.the-cryosphere.net/5/809/2011/) 453 ! equivalent to Arthern et al. 2010 eq. (4) "Nabarro-Herring" (doi:10.1029/2009JF001306) 454 ! Integration of the steady state equation 455 ! ln_smb approximated as a function of temperature 456 ln_smb = max((mean_temp(ikl) - TfSnow) * 5. / 60. + 8., 3.) 457 ! alpha0, alpha1 : correction coefficient as a function of ln_SMB from Ligtenberg 2011, adjusted for alpha1 458 alpha0 = max(1.435 - 0.151 * ln_smb, 0.25) 459 alpha1 = max(2.0111 - 0.2051 * ln_smb, 0.25) 460 E0 = C0 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha0 461 E1 = C1 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha1 462 z550 = log((rho_ice/mean_dens(ikl) - 1.)/(rho_ice/550. - 1.)) / E0 463 rho0 = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice 464 rho1 = exp(E1 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E1 * depth)) * rho_ice 465 if (depth <= z550) then 466 ro__SV(ikl, isl) = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice 467 else 468 ro__SV(ikl, isl) = exp(E1 * (depth - z550)) / (rho_ice / 550. - 1 + exp(E1 * (depth - z550))) * rho_ice 469 end if 470 depth = depth - dzsnSV(ikl,isl) / 2. 471 472 end do 473 474 END DO 475 476 END IF 477 478 479 ! + Numerics paramaters, SISVAT_ini 480 ! + ---------------------- 481 CALL SISVAT_ini(knon) 482 483 484 ! +--Read restart file 485 ! + ================================================= 486 487 INQUIRE(FILE = "startsis.nc", EXIST = file_exists) 488 IF (file_exists) THEN 489 CALL sisvatetat0("startsis.nc", ikl2i) 490 END IF 491 492 493 494 ! +--Output ascii file 495 ! + ================================================= 496 497 ! open output file 498 IF (ok_outfor) THEN 499 open(unit = un_outfor, status = 'replace', file = fn_outfor) 500 ikl = gp_outfor ! index sur la grille land ice 501 write(un_outfor, *) fn_outfor, ikl, dt__SV, rlon(ikl2i(ikl)), rlat(ikl2i(ikl)) 502 write(un_outfor, *) 'nsnow - albedo - z0m - z0h , dz [m,30], temp [K,41], rho [kg/m3,41], eta [kg/kg,41] & 503 & G1 [-,30], G2 [-,30], agesnow [d,30], history [-,30], DOP [m,30]' 504 END IF 505 506 END IF ! firstcall 507 ! + 508 ! + +++ INITIALISATION: END +++ 509 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 510 511 512 513 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 514 ! + READ FORCINGS 515 ! + ------------------------ 516 517 ! + Update Forcings for SISVAT given by the LMDZ model. 518 ! + 519 DO ikl = 1, knon 520 521 ! +--Atmospheric Forcing (INPUT) 522 ! + ^^^^^^^^^^^^^^^^^^^ ^^^^^ 523 za__SV(ikl) = zsl_height(ikl) ! surface layer height (fisr model level) [m] 524 Ua_min = 0.2 * sqrt(za__SV(ikl)) ! 525 VV__SV(ikl) = max(Ua_min, wind_velo(ikl)) ! Wind velocity [m/s] 526 TaT_SV(ikl) = temp_air(ikl) ! BL top Temperature [K] 527 ExnrSV(ikl) = pexner(ikl) ! Exner potential 528 rhT_SV(ikl) = dens_air(ikl) ! Air density 529 QaT_SV(ikl) = spechum(ikl) ! Specific humidity 530 ps__SV(ikl) = ps(ikl) ! surface pressure [Pa] 531 p1l_SV(ikl) = p1lay(ikl) ! lowest atm. layer press[Pa] 532 533 ! +--Surface properties 534 ! + ^^^^^^^^^^^^^^^^^^ 535 536 Z0m_SV(ikl) = z0m(ikl) ! Moment.Roughn.L. 537 Z0h_SV(ikl) = z0h(ikl) ! Moment.Roughn.L. 538 539 ! +--Energy Fluxes (INPUT) 540 ! + ^^^^^^^^^^^^^ ^^^^^ 541 coszSV(ikl) = max(czemin, rmu0(ikl)) ! cos(zenith.Dist.) 542 sol_SV(ikl) = swdown(ikl) ! downward Solar 543 IRd_SV(ikl) = lwdown(ikl) ! downward IR 544 rsolSV(ikl) = radsol(ikl) ! surface absorbed rad. 545 546 ! +--Water Fluxes (INPUT) 547 ! + ^^^^^^^^^^^^^ ^^^^^ 548 drr_SV(ikl) = precip_rain(ikl) ! Rain fall rate [kg/m2/s] 549 dsn_SV(ikl) = precip_snow(ikl) ! Snow fall rate [kg/m2/s] 550 551 ! #BS dbs_SV(ikl) = blowSN(i,j,n) 552 ! dbs_SV = Maximum potential erosion amount [kg/m2] 553 ! => Upper bound for eroded snow mass 554 ! uss_SV(ikl) = SLussl(i,j,n) ! u*qs* (only for Tv in sisvatesbl.f) 555 ! #BS if(dsn_SV(ikl)>eps12.and.erprev(i,j,n).gt.eps9) then 556 ! #BS dsnbSV(ikl) =1.0-min(qsHY(i,j,kB) !BS neglib. at kb ~100 magl) 557 ! #BS. /max(qshy(i,j,mz),eps9),unun) 558 ! #BS dsnbSV(ikl) = max(dsnbSV(ikl),erprev(i,j,n)/dsn_SV(ikl)) 559 ! #BS dsnbSV(ikl) = max(0.,min(1.,dsnbSV(ikl))) 560 ! #BS else 561 ! #BS dsnbSV(ikl) = 0. 562 ! #BS endif 563 ! dsnbSV is the drift fraction of deposited snow updated in sisvat.f 564 ! will be used for characterizing the Buffer Layer 565 ! (see update of Bros_N, G1same, G2same, zroOLD, zroNEW) 566 ! #BS if(n==1) qbs_HY(i,j) = dsnbSV(ikl) 567 qsnoSV(ikl) = snow_cont_air(ikl) 568 569 570 571 ! +--Soil/BL (INPUT) 572 ! + ^^^^^^^ ^^^^^ 573 alb0SV(ikl) = alb_soil(ikl) ! Soil background Albedo 574 AcoHSV(ikl) = AcoefH(ikl) 575 BcoHSV(ikl) = BcoefH(ikl) 576 AcoQSV(ikl) = AcoefQ(ikl) 577 BcoQSV(ikl) = BcoefQ(ikl) 578 cdH_SV(ikl) = min(cdragh(ikl),cdmax) 579 cdM_SV(ikl) = min(cdragm(ikl),cdmax) 580 rcdmSV(ikl) = sqrt(cdM_SV(ikl)) 581 Us_min = 0.01 582 us__SV(ikl) = max(Us_min, ustar(ikl)) 583 ram_sv(ikl) = 1. / (cdM_SV(ikl) * max(VV__SV(ikl), eps6)) 584 rah_sv(ikl) = 1. / (cdH_SV(ikl) * max(VV__SV(ikl), eps6)) 585 586 ! +--Energy Fluxes (INPUT/OUTPUT) 587 ! + ^^^^^^^^^^^^^ ^^^^^^^^^^^^ 588 !IF (.not.firstcall) THEN 589 Tsrfsv(ikl) = tsurf(ikl) !hj 12 03 2010 590 cld_SV(ikl) = cloudf(ikl) ! Cloudiness 591 !END IF 592 593 END DO 594 595 ! 596 ! + +++ READ FORCINGS: END +++ 597 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 598 599 600 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 601 ! +--SISVAT EXECUTION 602 ! + ---------------- 603 604 call INLANDSIS(SnoMod, BloMod, 1) 605 606 607 608 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 609 ! + RETURN RESULTS 610 ! + -------------- 611 ! + Return (compressed) SISVAT variables to LMDZ 612 ! + 613 DO ikl = 1, knon ! use only 1:knon (actual ice sheet..) 614 dflux_s(ikl) = dSdTSV(ikl) ! Sens.H.Flux T-Der. 615 dflux_l(ikl) = dLdTSV(ikl) ! Latn.H.Flux T-Der. 616 fluxsens(ikl) = HSs_sv(ikl) ! HS 617 fluxlat(ikl) = HLs_sv(ikl) ! HL 618 evap(ikl) = -1*HLs_sv(ikl) / LHvH2O ! Evaporation 619 erod(ikl) = 0. 620 621 IF (BloMod) THEN 622 ! + Blowing snow 623 624 ! SLussl(i,j,n)= 0. 625 ! #BS SLussl(i,j,n)= !Effective erosion 626 ! #BS. (- dbs_ER(ikl))/(dt*rhT_SV(ikl)) !~u*qs* from previous time step 627 ! #BS blowSN(i,j,n)= dt*uss_SV(ikl) !New max. pot. Erosion [kg/m2] 628 ! #BS. *rhT_SV(ikl) !(further bounded in sisvat_bsn.f) 629 ! #BS erprev(i,j,n) = dbs_Er(ikl)/dt__SV 630 erod(ikl) = dbs_Er(ikl) / dt__SV 631 ENDIF 632 633 ! + Check snow thickness, substract if too thick, add if too thin 634 635 sissnow(ikl) = 0. !() 636 DO isn = 1, isnoSV(ikl) 637 sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn) 638 END DO 639 640 IF (sissnow(ikl) .LE. sn_low) THEN !add snow 641 IF (isnoSV(ikl).GE.1) THEN 642 dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi) 643 toicSV(ikl) = toicSV(ikl) - sn_add 644 ELSE 645 write(*, *) 'Attention, bare ice... point ', ikl 646 isnoSV(ikl) = 1 647 istoSV(ikl, 1) = 0 648 ro__SV(ikl, 1) = 350. 649 dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi) ! 1. 650 eta_SV(ikl, 1) = epsi 651 TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2) 652 G1snSV(ikl, 1) = 0. 653 G2snSV(ikl, 1) = 0.3 654 agsnSV(ikl, 1) = 10. 655 toicSV(ikl) = toicSV(ikl) - sn_add 656 END IF 657 END IF 658 659 IF (sissnow(ikl) .ge. sn_upp) THEN !thinnen snow layer below 660 dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div 661 toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div 662 END IF 663 664 sissnow(ikl) = 0. 665 qsnow(ikl) = 0. 666 snow(ikl) = 0. 667 snowhgt(ikl) = 0. 668 669 DO isn = 1, isnoSV(ikl) 670 sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn) 671 snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn) 672 ! Etienne: check calc qsnow 673 qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn) 674 END DO 675 676 zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0) 677 ! Etienne: comment following line 678 ! snow(ikl) = sissnow(ikl)+toicSV(ikl) 679 snow(ikl) = sissnow(ikl) 680 681 to_ice(ikl) = toicSV(ikl) 682 runoff_lic(ikl) = RnofSV(ikl) ! RunOFF: intensity (flux due to melting + liquid precip) 683 fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing 684 ffonte(ikl)=fqfonte(ikl)*Lf_H2O 685 686 qsol(ikl) = 0. 687 DO isl = -nsol, 0 688 tsoil(ikl, 1 - isl) = TsisSV(ikl, isl) ! Soil Temperature 689 ! Etienne: check calc qsol 690 qsol(ikl) = qsol(ikl) & 691 + eta_SV(ikl, isl) * dz_dSV(isl) 692 END DO 693 agesno(ikl) = agsnSV(ikl, isnoSV(ikl)) ! [day] 694 695 alb1(ikl) = alb1sv(ikl) ! Albedo VIS 696 ! alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl) & 697 ! & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1 698 alb2(ikl)=alb2sv(ikl) 699 ! Albedo NIR 700 alb3(ikl) = alb3sv(ikl) ! Albedo FIR 701 ! 6 band Albedo 702 alb6(ikl,:)=alb6sv(ikl,:) 703 704 tsurf_new(ikl) = Tsrfsv(ikl) 705 706 qsurf(ikl) = QsT_SV(ikl) 707 emis_new(ikl) = eps0SL(ikl) 708 z0m(ikl) = Z0m_SV(ikl) 709 z0h(ikl) = Z0h_SV(ikl) 710 247 711 248 712 END DO 249 713 250 251 252 253 DO ikl=1,knon 254 255 256 ! Initialise variables 257 258 ispiSV(ikl) = 0 259 iiceSV(ikl) = 0 260 rusnSV(ikl) = 0. 261 toicSV(ikl) = 0. 262 isnoSV(ikl) = 0. ! # snow layers 263 istoSV(ikl,:) = 0. 264 eta_SV(ikl,:) = 0. 265 TsisSV(ikl,:) = 0. 266 ro__SV(ikl,:) = 0. 267 G1snSV(ikl,:) = 0. 268 G2snSV(ikl,:) = 0. 269 agsnSV(ikl,:) = 0. 270 dzsnSV(ikl,:) = 0. 271 zzsnsv(ikl,:) = 0. 272 BufsSV(ikl) = 0. 273 qsnoSV(ikl) = 0. ! BL snow content 274 zWEcSV(ikl) = 0. 275 dbs_SV(ikl) = 0. 276 dsnbSV(ikl) = 0. 277 esnbSV(ikl) = 0. 278 BrosSV(ikl) = 0. 279 BG1sSV(ikl) = 0. 280 BG2sSV(ikl) = 0. 281 SWS_SV(ikl) = 0. 282 RnofSV(ikl) = 0. ! RunOFF Intensity 283 RRs_SV(ikl) = 0. 284 DDs_SV(ikl) = 0. 285 VVs_SV(ikl) = 0. 286 cld_SV(ikl) = 0. 287 uts_SV(ikl) = 0. ! u*T* arbitrary 288 uqs_SV(ikl) = 0. ! u*q* " 289 uss_SV(ikl) = 0. ! u*s* " 290 LMO_SV(ikl) = 0. 291 292 293 ! Set variables 294 295 LSmask(ikl) = 1 ! Land/Sea Mask 296 isotSV(ikl) = 12 ! Soil Type -> 12= ice 297 iWaFSV(ikl) = 1 ! Soil Drainage 298 eps0SL(ikl )= 1. 299 alb0SV(ikl) = alb_soil(ikl) ! Soil Albedo 300 Z0m_SV(ikl) = z0m(ikl) ! Moment.Roughn.L. 301 Z0h_SV(ikl) = z0h(ikl) ! heat Roughn.L. 302 303 ! + Soil Upward IR Flux, Water Fluxes, roughness length 304 IRs_SV(ikl) = & 305 -eps0SL(ikl)* StefBo*(temp_air(ikl)**4) ! Upward IR Flux 306 Tsf_SV(ikl) = min(temp_air(ikl),TfSnow) 307 308 ! + Soil 309 DO isl = -nsol,0 310 TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2) !temp_air(ikl) !tsoil(ikl,1-isl) Soil Temperature 311 !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2) 312 eta_SV(ikl,isl) = epsi !etasoil(ikl,1-isl) Soil Water[m3/m3] 313 ro__SV(ikl,isl) = rhoIce !rosoil(ikl,1-isl) volumic mass 314 END DO 315 316 317 318 !! Initialise with snow 319 ! G1snSV(ikl,0) = 0. ! [-] 320 ! G2snSV(ikl,0) = 1.6 ! [-] [0.0001 m] 321 ! dzsnSV(ikl,0) = dz_dSV(0) ! [m] 322 323 324 ! if (snow(ikl) .GT. 0.) then 325 ! isnoSV(ikl) = 1 ! snow layers 326 ! istoSV(ikl,1:nsno) = 0 ! 0,...,5 : Snow History (see istdSV data) 327 ! eta_SV(ikl,1:nsno) = epsi 328 ! TsisSV(ikl,1:nsno) = tsoil(ikl,1) 329 ! ro__SV(ikl,1:nsno) = 350.0 330 ! G1snSV(ikl,1:nsno) = 0. ! [-] 331 ! G2snSV(ikl,1:nsno) = 1.6 ! [-] [0.0001 m] 332 ! agsnSV(ikl,1:nsno) = 50. ! [day] 333 ! dzsnSV(ikl,1) = snow(ikl)/max(ro__SV(ikl,1),epsi) ![m] 334 ! ! ecrete si trop de neige: 335 ! IF (snow(ikl) .ge. sn_upp) THEN !thinnen snow layer below 336 ! dzsnSV(ikl,1) = dzsnSV(ikl,1)/sn_div 337 ! toicSV(ikl) = toicSV(ikl)+dzsnSV(ikl,1)*ro__SV(ikl,1)/sn_div 338 ! END IF 339 ! zzsnsv(ikl,1) = dzsnSV(ikl,1) ! Total snow pack thickness 340 ! endif 341 342 343 ! Initialise la neige avec un profil de densité prochde des conditions de Dôme C (~10m de neige avec 19 niveaux) (Etienne): 344 isnoSV(ikl) = 19 345 istoSV(ikl,1:isnoSV(ikl)) = 100 346 ro__SV(ikl,1:isnoSV(ikl)) = 350. 347 eta_SV(ikl,1:isnoSV(ikl)) = epsi 348 TsisSV(ikl,1:isnoSV(ikl)) = min(tsoil(ikl,1),TfSnow-0.2) 349 G1snSV(ikl,1:isnoSV(ikl)) = 0 350 G2snSV(ikl,1:isnoSV(ikl)) = 1.6 351 agsnSV(ikl,1:isnoSV(ikl)) = 50. 352 dzsnSV(ikl,19) = 0.015 353 dzsnSV(ikl,18) =0.015 354 dzsnSV(ikl,17) =0.020 355 dzsnSV(ikl,16) =0.030 356 dzsnSV(ikl,15) =0.040 357 dzsnSV(ikl,14) =0.060 358 dzsnSV(ikl,13) =0.080 359 dzsnSV(ikl,12) =0.110 360 dzsnSV(ikl,11) =0.150 361 dzsnSV(ikl,10) =0.200 362 dzsnSV(ikl,9) =0.300 363 dzsnSV(ikl,8) =0.420 364 dzsnSV(ikl,7) =0.780 365 dzsnSV(ikl,6) =1.020 366 dzsnSV(ikl,5) =0.980 367 dzsnSV(ikl,4) =1.020 368 dzsnSV(ikl,3) =3.970 369 dzsnSV(ikl,2) =1.020 370 dzsnSV(ikl,1) =0.100 371 372 373 END DO 374 375 ! +--Surface Fall Line Slope 376 ! + ----------------------- 377 IF (SnoMod) THEN 378 DO ikl=1,knon 379 slopSV(ikl) = slope(ikl) 380 SWf_SV(ikl) = & ! Normalized Decay of the 381 exp(-dt__SV & ! Surficial Water Content 382 /(c1_zuo & !(Zuo and Oerlemans 1996, 383 +c2_zuo*exp(-c3_zuo*abs(slopSV(ikl))))) ! J.Glacio. 42, 305--317) 384 END DO 385 END IF 386 387 ! + SISVAT_ini (as for use with MAR, but not computing soil layers) 388 ! + ------------------------------------------------------------- 389 ! write(*,'(/a)') 'Start SISVAT initialization: SISVAT_ini' 390 CALL SISVAT_ini(knon) 391 392 393 ! +--Read restart file 394 ! + ================================================= 395 396 INQUIRE(FILE="startsis.nc", EXIST=file_exists) 397 IF (file_exists) THEN 398 CALL sisvatetat0("startsis.nc",ikl2i) 714 IF (ok_outfor) THEN 715 ikl= gp_outfor 716 write(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++' 717 write(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl) 718 write(un_outfor, *) dzsnSV(ikl, :) 719 write(un_outfor, *) TsisSV(ikl, :) 720 write(un_outfor, *) ro__SV(ikl, :) 721 write(un_outfor, *) eta_SV(ikl, :) 722 write(un_outfor, *) G1snSV(ikl, :) 723 write(un_outfor, *) G2snSV(ikl, :) 724 write(un_outfor, *) agsnSV(ikl, :) 725 write(un_outfor, *) istoSV(ikl, :) 726 write(un_outfor, *) DOPsnSV(ikl, :) 727 ENDIF 728 729 730 731 ! + ----------------------------- 732 ! + END --- RETURN RESULTS 733 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 734 IF (lafin) THEN 735 fichnom = "restartsis.nc" 736 CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat) 737 738 IF (ok_outfor) THEN 739 close(unit = un_outfor) 740 END IF 399 741 END IF 400 401 402 403 ! +--Output ascii file 404 ! + ================================================= 405 406 407 408 ! open output file 409 IF (ok_outfor) THEN 410 open(unit=un_outfor,status='replace',file=fn_outfor) 411 ikl=gp_outfor ! index sur la grille land ice 412 write(un_outfor,*) fn_outfor, ikl, dt__SV 413 write(un_outfor,*) 'nsnow - albedo - z0m - z0h , dz [m,35], temp [K,46], rho [kg/m3,46], eta [kg/kg,46] & 414 & G1 [-,35], G2 [-,35], agesnow [d,35], history [-,35]' 415 416 END IF 417 418 END IF ! firstcall 419 ! + 420 ! + +++ INITIALISATION: END +++ 421 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 422 423 424 425 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 426 ! + READ FORCINGS 427 ! + ------------------------ 428 429 ! + Update Forcings for SISVAT given by the LMDZ model. 430 ! + 431 DO ikl=1,knon 432 433 ! +--Atmospheric Forcing (INPUT) 434 ! + ^^^^^^^^^^^^^^^^^^^ ^^^^^ 435 zSBLSV = 1000. ! [m] 436 za__SV(ikl) = zsl_height(ikl) ! surface layer height (fisr model level) [m] 437 Ua_min = epsi ! 438 Ua_min = 0.2 * sqrt(za__SV(ikl) ) ! 439 VV__SV(ikl) = max(Ua_min, wind_velo(ikl)) ! Wind velocity [m/s] 440 TaT_SV(ikl) = temp_air(ikl) ! BL top Temperature [K] 441 ExnrSV(ikl) = pexner(ikl) ! Exner potential 442 rhT_SV(ikl) = dens_air(ikl) ! Air density 443 QaT_SV(ikl) = spechum(ikl) ! Specific humidity 444 ps__SV(ikl) = ps(ikl) ! surface pressure [Pa] 445 p1l_SV(ikl) = p1lay(ikl) ! lowest atm. layer press[Pa] 446 447 ! +--Surface properties 448 ! + ^^^^^^^^^^^^^^^^^^ 449 450 Z0m_SV(ikl) = z0m(ikl) ! Moment.Roughn.L. 451 Z0h_SV(ikl) = z0h(ikl) ! Moment.Roughn.L. 452 453 ! +--Energy Fluxes (INPUT) 454 ! + ^^^^^^^^^^^^^ ^^^^^ 455 coszSV(ikl) = max(czemin,rmu0(ikl)) ! cos(zenith.Dist.) 456 sol_SV(ikl) = swdown(ikl) ! downward Solar 457 IRd_SV(ikl) = lwdown(ikl) ! downward IR 458 rsolSV(ikl) = radsol(ikl) ! surface absorbed rad. 459 460 ! +--Water Fluxes (INPUT) 461 ! + ^^^^^^^^^^^^^ ^^^^^ 462 drr_SV(ikl) = precip_rain(ikl) ! Rain fall rate [kg/m2/s] 463 dsn_SV(ikl) = precip_snow(ikl) ! Snow fall rate [kg/m2/s] 464 !c #BS dbsnow = -SLussl(i,j,n) ! Erosion 465 !c #BS. *dtPhys *rhT_SV(ikl) /ro_Wat 466 !c #BS dsnbSV(ikl) = snow_adv(ikl) ! min(max(zero,dbsnow) 467 !c #BS. / max(epsi,d_snow),unun) 468 !c #BS dbs_SV(ikl) = snow_cont_air(ikl) 469 !c #BS blowSN(i,j,n) ! [kg/m2] 470 471 ! +--Soil/BL (INPUT) 472 ! + ^^^^^^^ ^^^^^ 473 alb0SV(ikl) = alb_soil(ikl) ! Soil background Albedo 474 AcoHSV(ikl) = AcoefH(ikl) 475 BcoHSV(ikl) = BcoefH(ikl) 476 AcoQSV(ikl) = AcoefQ(ikl) 477 BcoQSV(ikl) = BcoefQ(ikl) 478 cdH_SV(ikl) = cdragh(ikl) 479 cdM_SV(ikl) = cdragm(ikl) 480 Us_min = 0.01 481 us__SV(ikl) = max(Us_min, ustar(ikl)) 482 ram_sv(ikl) = 1./(cdragm(ikl)*max(VV__SV(ikl),eps6)) 483 rah_sv(ikl) = 1./(cdragh(ikl)*max(VV__SV(ikl),eps6)) 484 485 ! +--Energy Fluxes (INPUT/OUTPUT) 486 ! + ^^^^^^^^^^^^^ ^^^^^^^^^^^^ 487 IF (.not.firstcall) THEN 488 Tsf_SV(ikl) = tsurf(ikl) !hj 12 03 2010 489 cld_SV(ikl) = cloudf(ikl) ! Cloudiness 490 END IF 491 492 493 END DO 494 495 ! 496 ! + +++ READ FORCINGS: END +++ 497 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 498 499 500 501 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 502 ! +--SISVAT EXECUTION 503 ! + ---------------- 504 505 call INLANDSIS(SnoMod,BloMod,1) 506 507 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 508 ! + RETURN RESULTS 509 ! + -------------- 510 ! + Return (compressed) SISVAT variables to LMDZ 511 ! + 512 DO ikl=1,knon ! use only 1:knon (actual ice sheet..) 513 runoff_lic(ikl) = RnofSV(ikl)*dtime ! RunOFF: intensity* time step 514 dflux_s(ikl) = dSdTSV(ikl) ! Sens.H.Flux T-Der. 515 dflux_l(ikl) = dLdTSV(ikl) ! Latn.H.Flux T-Der. 516 fluxsens(ikl) = HSs_sv(ikl) ! HS 517 fluxlat(ikl) = HLs_sv(ikl) ! HL 518 evap(ikl) = HLs_sv(ikl)/LHvH2O ! Evaporation 519 snow(ikl) = 0. 520 snowhgt(ikl) = 0. 521 qsnow(ikl) = 0. 522 qsol(ikl) = 0. 523 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 524 ! + 525 ! + Check snow thickness, substract if too thick (commended by etienne: add if too thin) 526 527 sissnow(ikl) = 0. !() 528 DO isn = 1,isnoSV(ikl) 529 sissnow(ikl) = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn) 530 END DO 531 532 IF (sissnow(ikl) .LE. sn_low) THEN !add snow 533 IF (isnoSV(ikl).GE.1) THEN 534 dzsnSV(ikl,1) = dzsnSV(ikl,1) + sn_add/max(ro__SV(ikl,1),epsi) 535 toicSV(ikl) = toicSV(ikl) - sn_add 536 ! ELSE 537 ! write(*,*) 'Attention, bare ice... point ',ikl 538 ! isnoSV(ikl) = 1 539 ! istoSV(ikl,1) = 0 540 ! ro__SV(ikl,1) = 350. 541 ! dzsnSV(ikl,1) = sn_add/max(ro__SV(ikl,1),epsi) ! 1. 542 ! eta_SV(ikl,1) = epsi 543 ! TsisSV(ikl,1) = min(TsisSV(ikl,0),TfSnow-0.2) 544 ! G1snSV(ikl,1) = 0. 545 ! G2snSV(ikl,1) = 0.3 546 ! agsnSV(ikl,1) = 10. 547 ! toicSV(ikl) = toicSV(ikl) - sn_add 548 END IF 549 END IF 550 551 IF (sissnow(ikl) .ge. sn_upp) THEN !thinnen snow layer below 552 dzsnSV(ikl,1) = dzsnSV(ikl,1)/sn_div 553 toicSV(ikl) = toicSV(ikl)+dzsnSV(ikl,1)*ro__SV(ikl,1)/sn_div 554 END IF 555 556 sissnow(ikl) = 0. !() 557 558 DO isn = 1,isnoSV(ikl) 559 sissnow(ikl) = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn) 560 snowhgt(ikl) = snowhgt(ikl)+dzsnSV(ikl,isn) 561 qsnow(ikl) = qsnow(ikl)+1e03*eta_SV(ikl,isn)*dzsnSV(ikl,isn) 562 END DO 563 564 ! Etienne: pourquoi ajouter toicSV ici? Pour bilan d'eau? 565 snow(ikl) = sissnow(ikl)+toicSV(ikl) 566 to_ice(ikl) = toicSV(ikl) 567 568 569 DO isl = -nsol,0 570 tsoil(ikl,1-isl) = TsisSV(ikl,isl) ! Soil Temperature 571 qsol(ikl) = qsol(ikl) & 572 +eta_SV(ikl,isl) * dz_dSV(isl) 573 END DO 574 agesno(ikl) = agsnSV(ikl,isnoSV(ikl)) ! [day] 575 576 alb1(ikl) = alb1sv(ikl) ! Albedo VIS 577 alb2(ikl) = ((So1dSV-f1)*alb1sv(ikl) & 578 & +So2dSV*alb2sv(ikl)+So3dSV*alb3sv(ikl))/f1 579 ! Albedo NIR 580 alb3(ikl) = alb3sv(ikl) ! Albedo FIR 581 582 tsurf_new(ikl) =Tsrfsv(ikl) 583 584 zfra(ikl) = max(min(isnoSV(ikl)-iiceSV(ikl),1),0) 585 qsurf(ikl) = QaT_SV(ikl) 586 emis_new(ikl) = eps0SL(ikl) 587 z0m(ikl) = Z0m_SV(ikl) 588 z0h(ikl) = Z0h_SV(ikl) 589 590 END DO ! ikl 591 592 593 594 595 596 597 ! write variables in output file 598 599 IF (ok_outfor) THEN 600 ikl=gp_outfor 601 602 ! write(un_outfor,*) 'nsnow [-,1], dz [m,35], temp [K,46], rho [kg/m3,46], eta [kg/kg,46]' 603 ! write(un_outfor,*) 'G1 [-,35], G2 [-,35], agesnow [d,35], history [-,35]' 604 write(un_outfor,*) '+++++++++++++++++++++++++++++++++++++++++++++++' 605 write(un_outfor,*) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl) 606 write(un_outfor,*) dzsnSV(ikl,:) 607 write(un_outfor,*) TsisSV(ikl,:) 608 write(un_outfor,*) ro__SV(ikl,:) 609 write(un_outfor,*) eta_SV(ikl,:) 610 write(un_outfor,*) G1snSV(ikl,:) 611 write(un_outfor,*) G2snSV(ikl,:) 612 write(un_outfor,*) agsnSV(ikl,:) 613 write(un_outfor,*) istoSV(ikl,:) 614 615 ENDIF 616 617 618 619 620 ! + ----------------------------- 621 ! + END --- RETURN RESULTS 622 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 623 IF (lafin) THEN 624 fichnom = "restartsis.nc" 625 CALL sisvatredem("restartsis.nc",ikl2i,rlon,rlat) 626 627 IF (ok_outfor) THEN 628 close(unit=un_outfor) 629 END IF 630 END IF 631 632 633 END SUBROUTINE surf_inlandsis 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 !======================================================================= 649 650 SUBROUTINE get_soil_levels(dz1, dz2, lambda) 651 ! ====================================================================== 652 ! Routine to compute the vertical discretization of the soil in analogy 653 ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case 654 ! of SISVAT, therefore it's needed here. 655 ! 656 USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root 657 USE mod_phys_lmdz_para 658 USE VAR_SV 659 660 661 ! INCLUDE "dimsoil.h" 662 663 REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1 664 REAL, INTENT(OUT) :: lambda 665 666 667 !----------------------------------------------------------------------- 668 ! Depthts: 669 ! -------- 670 REAL fz,rk,fz1,rk1,rk2 671 REAL min_period, dalph_soil 672 INTEGER ierr,jk 673 674 fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) 675 676 ! write(*,*)'Start soil level computation' 677 !----------------------------------------------------------------------- 678 ! Calculation of some constants 679 ! NB! These constants do not depend on the sub-surfaces 680 !----------------------------------------------------------------------- 681 !----------------------------------------------------------------------- 682 ! ground levels 683 ! grnd=z/l where l is the skin depth of the diurnal cycle: 684 !----------------------------------------------------------------------- 685 686 min_period=1800. ! en secondes 687 dalph_soil=2. ! rapport entre les epaisseurs de 2 couches succ. 688 ! !$OMP MASTER 689 ! IF (is_mpi_root) THEN 690 ! OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr) 691 ! IF (ierr == 0) THEN ! Read file only if it exists 692 ! READ(99,*) min_period 693 ! READ(99,*) dalph_soil 694 ! PRINT*,'Discretization for the soil model' 695 ! PRINT*,'First level e-folding depth',min_period, & 696 ! ' dalph',dalph_soil 697 ! CLOSE(99) 698 ! END IF 699 ! ENDIF 700 ! !$OMP END MASTER 701 ! CALL bcast(min_period) 702 ! CALL bcast(dalph_soil) 703 704 ! la premiere couche represente un dixieme de cycle diurne 705 fz1=SQRT(min_period/3.14) 706 707 DO jk=1,nsoilmx 708 rk1=jk 709 rk2=jk-1 710 dz2(jk)=fz(rk1)-fz(rk2) 711 ENDDO 712 DO jk=1,nsoilmx-1 713 rk1=jk+.5 714 rk2=jk-.5 715 dz1(jk)=1./(fz(rk1)-fz(rk2)) 716 ENDDO 717 lambda=fz(.5)*dz1(1) 718 PRINT*,'full layers, intermediate layers (seconds)' 719 DO jk=1,nsoilmx 720 rk=jk 721 rk1=jk+.5 722 rk2=jk-.5 723 PRINT *,'fz=', & 724 fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14 725 ENDDO 726 727 END SUBROUTINE get_soil_levels 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 !=========================================================================== 746 747 SUBROUTINE SISVAT_ini(knon) 748 749 !C +------------------------------------------------------------------------+ 750 !C | MAR SISVAT_ini Jd 11-10-2007 MAR | 751 !C | SubRoutine SISVAT_ini generates non time dependant SISVAT parameters | 752 !C +------------------------------------------------------------------------+ 753 !C | PARAMETERS: klonv: Total Number of columns = | 754 !C | ^^^^^^^^^^ = Total Number of continental grid boxes | 755 !C | X Number of Mosaic Cell per grid box | 756 !C | | 757 !C | INPUT: dt__SV : Time Step [s] | 758 !C | ^^^^^ dz_dSV : Layer Thickness [m] | 759 !C | | 760 !C | OUTPUT: [-] | 761 !C | ^^^^^^ rocsSV : Soil Contrib. to (ro c)_s exclud.Water [J/kg/K] | 762 !C | etamSV : Soil Minimum Humidity [m3/m3] | 763 !C | (based on a prescribed Soil Relative Humidity) | 764 !C | s1__SV : Factor of eta**( b+2) in Hydraul.Diffusiv. | 765 !C | s2__SV : Factor of eta**( b+2) in Hydraul.Conduct. | 766 !C | aKdtSV : KHyd: Piecewise Linear Profile: a * dt [m] | 767 !C | bKdtSV : KHyd: Piecewise Linear Profile: b * dt [m/s] | 768 !C | dzsnSV(0): Soil first Layer Thickness [m] | 769 !C | dzmiSV : Distance between two contiguous levels [m] | 770 !C | dz78SV : 7/8 (Layer Thickness) [m] | 771 !C | dz34SV : 3/4 (Layer Thickness) [m] | 772 !C | dz_8SV : 1/8 (Layer Thickness) [m] | 773 !C | dzAvSV : 1/8 dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1) [m] | 774 !C | dtz_SV : dt/dz [s/m] | 775 !C | OcndSV : Swab Ocean / Soil Ratio [-] | 776 !C | Implic : Implicit Parameter (0.5: Crank-Nicholson) | 777 !C | Explic : Explicit Parameter = 1.0 - Implic | 778 !C | | 779 !C | # OPTIONS: #ER: Richards Equation is not smoothed | 780 !C | # ^^^^^^^ #kd: De Ridder Discretization | 781 !C | # #SH: Hapex-Sahel Values ! 782 !C | | 783 !C +------------------------------------------------------------------------+ 784 ! 785 ! 786 787 !C +--Global Variables 788 !C + ================ 789 790 USE dimphy 791 USE VARphy 792 USE VAR_SV 793 USE VARdSV 794 USE VAR0SV 795 USE VARxSV 796 USE VARtSV 797 USE VARxSV 798 USE VARySV 799 IMPLICIT NONE 800 801 802 803 !C +--Arguments 804 !C + ================== 805 INTEGER,INTENT(IN) :: knon 806 807 !C +--Internal Variables 808 !C + ================== 809 810 INTEGER :: ivt ,ist ,ikl ,isl ,isn ,ikh 811 INTEGER :: misl_2,nisl_2 812 REAL :: d__eta,eta__1,eta__2,Khyd_1,Khyd_2 813 REAL,PARAMETER :: RHsMin= 0.001 ! Min.Soil Relative Humidity 814 REAL :: PsiMax ! Max.Soil Water Potential 815 REAL :: a_Khyd,b_Khyd ! Piecewis.https://www.lequipe.fr/Water Conductivity 816 817 818 !c #WR REAL :: Khyd_x,Khyd_y 819 820 821 822 !C +--Non Time Dependant SISVAT parameters 823 !C + ==================================== 824 825 !C +--Soil Discretization 826 !C + ------------------- 827 828 !C +--Numerical Scheme Parameters 829 !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^ 742 743 END SUBROUTINE surf_inlandsis 744 745 746 !======================================================================= 747 748 SUBROUTINE get_soil_levels(dz1, dz2, lambda) 749 ! ====================================================================== 750 ! Routine to compute the vertical discretization of the soil in analogy 751 ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case 752 ! of SISVAT, therefore it's needed here. 753 ! 754 USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root 755 USE mod_phys_lmdz_para 756 USE VAR_SV 757 758 759 ! INCLUDE "dimsoil.h" 760 761 REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1 762 REAL, INTENT(OUT) :: lambda 763 764 765 !----------------------------------------------------------------------- 766 ! Depthts: 767 ! -------- 768 REAL fz, rk, fz1, rk1, rk2 769 REAL min_period, dalph_soil 770 INTEGER ierr, jk 771 772 fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.) 773 774 ! write(*,*)'Start soil level computation' 775 !----------------------------------------------------------------------- 776 ! Calculation of some constants 777 ! NB! These constants do not depend on the sub-surfaces 778 !----------------------------------------------------------------------- 779 !----------------------------------------------------------------------- 780 ! ground levels 781 ! grnd=z/l where l is the skin depth of the diurnal cycle: 782 !----------------------------------------------------------------------- 783 784 min_period = 1800. ! en secondes 785 dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ. 786 ! !$OMP MASTER 787 ! IF (is_mpi_root) THEN 788 ! OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr) 789 ! IF (ierr == 0) THEN ! Read file only if it exists 790 ! READ(99,*) min_period 791 ! READ(99,*) dalph_soil 792 ! PRINT*,'Discretization for the soil model' 793 ! PRINT*,'First level e-folding depth',min_period, & 794 ! ' dalph',dalph_soil 795 ! CLOSE(99) 796 ! END IF 797 ! ENDIF 798 ! !$OMP END MASTER 799 ! CALL bcast(min_period) 800 ! CALL bcast(dalph_soil) 801 802 ! la premiere couche represente un dixieme de cycle diurne 803 fz1 = SQRT(min_period / 3.14) 804 805 DO jk = 1, nsoilmx 806 rk1 = jk 807 rk2 = jk - 1 808 dz2(jk) = fz(rk1) - fz(rk2) 809 ENDDO 810 DO jk = 1, nsoilmx - 1 811 rk1 = jk + .5 812 rk2 = jk - .5 813 dz1(jk) = 1. / (fz(rk1) - fz(rk2)) 814 ENDDO 815 lambda = fz(.5) * dz1(1) 816 DO jk = 1, nsoilmx 817 rk = jk 818 rk1 = jk + .5 819 rk2 = jk - .5 820 ENDDO 821 822 END SUBROUTINE get_soil_levels 823 824 825 !=========================================================================== 826 827 SUBROUTINE SISVAT_ini(knon) 828 829 !C +------------------------------------------------------------------------+ 830 !C | MAR SISVAT_ini Jd 11-10-2007 MAR | 831 !C | SubRoutine SISVAT_ini generates non time dependant SISVAT parameters | 832 !C +------------------------------------------------------------------------+ 833 !C | PARAMETERS: klonv: Total Number of columns = | 834 !C | ^^^^^^^^^^ = Total Number of continental grid boxes | 835 !C | X Number of Mosaic Cell per grid box | 836 !C | | 837 !C | INPUT: dt__SV : Time Step [s] | 838 !C | ^^^^^ dz_dSV : Layer Thickness [m] | 839 !C | | 840 !C | OUTPUT: [-] | 841 !C | ^^^^^^ rocsSV : Soil Contrib. to (ro c)_s exclud.Water [J/kg/K] | 842 !C | etamSV : Soil Minimum Humidity [m3/m3] | 843 !C | (based on a prescribed Soil Relative Humidity) | 844 !C | s1__SV : Factor of eta**( b+2) in Hydraul.Diffusiv. | 845 !C | s2__SV : Factor of eta**( b+2) in Hydraul.Conduct. | 846 !C | aKdtSV : KHyd: Piecewise Linear Profile: a * dt [m] | 847 !C | bKdtSV : KHyd: Piecewise Linear Profile: b * dt [m/s] | 848 !C | dzsnSV(0): Soil first Layer Thickness [m] | 849 !C | dzmiSV : Distance between two contiguous levels [m] | 850 !C | dz78SV : 7/8 (Layer Thickness) [m] | 851 !C | dz34SV : 3/4 (Layer Thickness) [m] | 852 !C | dz_8SV : 1/8 (Layer Thickness) [m] | 853 !C | dzAvSV : 1/8 dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1) [m] | 854 !C | dtz_SV : dt/dz [s/m] | 855 !C | OcndSV : Swab Ocean / Soil Ratio [-] | 856 !C | Implic : Implicit Parameter (0.5: Crank-Nicholson) | 857 !C | Explic : Explicit Parameter = 1.0 - Implic | 858 !C | | 859 !C | # OPTIONS: #ER: Richards Equation is not smoothed | 860 !C | # ^^^^^^^ #kd: De Ridder Discretization | 861 !C | # #SH: Hapex-Sahel Values ! 862 !C | | 863 !C +------------------------------------------------------------------------+ 864 ! 865 ! 866 867 !C +--Global Variables 868 !C + ================ 869 870 USE dimphy 871 USE VARphy 872 USE VAR_SV 873 USE VARdSV 874 USE VAR0SV 875 USE VARxSV 876 USE VARtSV 877 USE VARxSV 878 USE VARySV 879 IMPLICIT NONE 880 881 882 883 !C +--Arguments 884 !C + ================== 885 INTEGER, INTENT(IN) :: knon 886 887 !C +--Internal Variables 888 !C + ================== 889 890 INTEGER :: ivt, ist, ikl, isl, isn, ikh 891 INTEGER :: misl_2, nisl_2 892 REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2 893 REAL, PARAMETER :: RHsMin = 0.001 ! Min.Soil Relative Humidity 894 REAL :: PsiMax ! Max.Soil Water Potential 895 REAL :: a_Khyd, b_Khyd ! Water conductivity 896 897 898 !c #WR REAL :: Khyd_x,Khyd_y 899 900 901 902 !C +--Non Time Dependant SISVAT parameters 903 !C + ==================================== 904 905 !C +--Soil Discretization 906 !C + ------------------- 907 908 !C +--Numerical Scheme Parameters 909 !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^ 830 910 Implic = 0.75 ! 0.5 <==> Crank-Nicholson 831 911 Explic = 1.00 - Implic ! 832 833 !C +--Soil/Snow Layers Indices 834 !C + ^^^^^^^^^^^^^^^^^^^^^^^^ 835 DO isl=-nsol,0 836 islpSV(isl) = isl+1 837 islpSV(isl) = min( islpSV(isl),0) 838 islmSV(isl) = isl-1 839 islmSV(isl) = max(-nsol,islmSV(isl)) 840 END DO 841 842 DO isn=1,nsno 843 isnpSV(isn) = isn+1 844 isnpSV(isn) = min( isnpSV(isn),nsno) 845 END DO 846 847 !C +--Soil Layers Thicknesses 848 !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 849 ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels! 850 !c #kd IF (nsol.gt.4) THEN 851 !c #kd DO isl=-5,-nsol,-1 852 !c #kd dz_dSV(isl)= 1. 853 !c #kd END DO 854 !c #kd END IF 855 ! 856 ! IF (nsol.ne.4) THEN 857 ! DO isl= 0,-nsol,-1 858 ! misl_2 = -mod(isl,2) 859 ! nisl_2 = -isl/2 860 ! dz_dSV(isl)=(((1-misl_2) * 0.001 861 ! . + misl_2 * 0.003) * 10**(nisl_2)) * 4. 862 !C +... dz_dSV(0) = Hapex-Sahel Calibration: 4 mm 863 ! 864 !c +SH dz_dSV(isl)=(((1-misl_2) * 0.001 865 !c +SH. + misl_2 * 0.003) * 10**(nisl_2)) * 1. 866 ! 867 !c #05 dz_dSV(isl)=(((1-misl_2) * 0.001 868 !c #05. + misl_2 * 0.008) * 10**(nisl_2)) * 0.5 869 ! END DO 870 ! dz_dSV(0) = 0.001 871 ! dz_dSV(-1) = dz_dSV(-1) - dz_dSV(0) + 0.004 872 ! END IF 873 874 875 zz_dSV = 0. 876 DO isl=-nsol,0 877 dzmiSV(isl) = 0.500*(dz_dSV(isl) +dz_dSV(islmSV(isl))) 878 dziiSV(isl) = 0.500* dz_dSV(isl) /dzmiSV(isl) 879 dzi_SV(isl) = 0.500* dz_dSV(islmSV(isl))/dzmiSV(isl) 880 dtz_SV(isl) = dt__SV /dz_dSV(isl) 881 dtz_SV2(isl) = 1. /dz_dSV(isl) 882 dz78SV(isl) = 0.875* dz_dSV(isl) 883 dz34SV(isl) = 0.750* dz_dSV(isl) 884 dz_8SV(isl) = 0.125* dz_dSV(isl) 885 dzAvSV(isl) = 0.125* dz_dSV(islmSV(isl)) & 886 & + 0.750* dz_dSV(isl) & 887 & + 0.125* dz_dSV(islpSV(isl)) 888 zz_dSV = zz_dSV+dz_dSV(isl) 889 END DO 890 DO ikl=1,knon !v 891 dzsnSV(ikl,0) = dz_dSV(0) 892 END DO 893 894 !C +--Conversion to a 50 m Swab Ocean Discretization 895 !C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 896 OcndSV = 0. 897 DO isl=-nsol,0 898 OcndSV = OcndSV +dz_dSV(isl) 899 END DO 900 OcndSV = 50. /OcndSV 901 902 903 !C +--Secondary Soil Parameters 904 !C + ------------------------------- 905 906 DO ist=0,nsot 907 rocsSV(ist)=(1.0-etadSV(ist))*1.2E+6 ! Soil Contrib. to (ro c)_s 908 s1__SV(ist)= bCHdSV(ist) & ! Factor of (eta)**(b+2) 909 & *psidSV(ist) *Ks_dSV(ist) & ! in DR97, Eqn.(3.36) 910 & /(etadSV(ist)**( bCHdSV(ist)+3.)) ! 911 s2__SV(ist)= Ks_dSV(ist) & ! Factor of (eta)**(2b+3) 912 & /(etadSV(ist)**(2.*bCHdSV(ist)+3.)) ! in DR97, Eqn.(3.35) 913 914 !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity) 915 !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 916 Psimax = -(log(RHsMin))/7.2E-5 ! DR97, Eqn 3.15 Inversion 917 etamSV(ist) = etadSV(ist) & 918 & *(PsiMax/psidSV(ist))**(-min(10.,1./bCHdSV(ist))) 919 END DO 920 etamSV(12) = 0. 921 922 !C +--Piecewise Hydraulic Conductivity Profiles 923 !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 924 DO ist=0,nsot 925 926 927 d__eta = etadSV(ist)/nkhy 928 eta__1 = 0. 929 eta__2 = d__eta 930 DO ikh=0,nkhy 931 Khyd_1 = s2__SV(ist) & ! DR97, Eqn.(3.35) 932 & *(eta__1 **(2. *bCHdSV(ist)+3.)) ! 933 Khyd_2 = s2__SV(ist) &! 934 & *(eta__2 **(2. *bCHdSV(ist)+3.)) ! 935 936 a_Khyd = (Khyd_2-Khyd_1)/d__eta ! 937 b_Khyd = Khyd_1-a_Khyd *eta__1 ! 938 !c #WR Khyd_x = a_Khyd*eta__1 +b_Khyd ! 939 !c #WR Khyd_y = a_Khyd*eta__2 +b_Khyd ! 940 aKdtSV(ist,ikh) = a_Khyd * dt__SV ! 941 bKdtSV(ist,ikh) = b_Khyd * dt__SV ! 942 943 eta__1 = eta__1 + d__eta 944 eta__2 = eta__2 + d__eta 945 END DO 946 END DO 947 948 949 return 950 951 END SUBROUTINE SISVAT_ini 952 953 954 955 956 957 958 959 !*************************************************************************** 960 961 SUBROUTINE sisvatetat0 (fichnom,ikl2i) 962 963 USE dimphy 964 USE mod_grid_phy_lmdz 965 USE mod_phys_lmdz_para 966 967 USE iostart 968 USE VAR_SV 969 USE VARdSV 970 USE VARxSV 971 USE VARtSV 972 USE indice_sol_mod 973 974 IMPLICIT none 975 !====================================================================== 976 ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009 977 ! Objet: Lecture du fichier de conditions initiales pour SISVAT 978 !====================================================================== 979 include "netcdf.inc" 980 ! include "indicesol.h" 981 982 ! include "dimsoil.h" 983 include "clesphys.h" 984 include "thermcell.h" 985 include "compbl.h" 986 987 !====================================================================== 988 CHARACTER(LEN=*) :: fichnom 989 990 991 INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i 992 REAL, DIMENSION(klon) :: rlon 993 REAL, DIMENSION(klon) :: rlat 994 995 ! les variables globales ecrites dans le fichier restart 996 REAL, DIMENSION(klon) :: isno 997 REAL, DIMENSION(klon) :: ispi 998 REAL, DIMENSION(klon) :: iice 999 REAL, DIMENSION(klon) :: rusn 1000 REAL, DIMENSION(klon, nsno) :: isto 1001 1002 REAL, DIMENSION(klon, nsismx) :: Tsis 1003 REAL, DIMENSION(klon, nsismx) :: eta 1004 REAL, DIMENSION(klon, nsismx) :: ro 1005 1006 REAL, DIMENSION(klon, nsno) :: dzsn 1007 REAL, DIMENSION(klon, nsno) :: G1sn 1008 REAL, DIMENSION(klon, nsno) :: G2sn 1009 REAL, DIMENSION(klon, nsno) :: agsn 1010 1011 REAL, DIMENSION(klon) :: toic 1012 1013 1014 INTEGER :: isl, ikl, i, isn , errT, erreta, errro, errdz, snopts 1015 CHARACTER (len=2) :: str2 1016 LOGICAL :: found 1017 1018 errT=0 1019 errro=0 1020 erreta=0 1021 errdz=0 1022 snopts=0 1023 ! Ouvrir le fichier contenant l'etat initial: 1024 1025 CALL open_startphy(fichnom) 1026 1027 ! Lecture des latitudes, longitudes (coordonnees): 1028 1029 CALL get_field("latitude",rlat,found) 1030 CALL get_field("longitude",rlon,found) 1031 1032 CALL get_field("n_snows", isno,found) 1033 IF (.NOT. found) THEN 1034 PRINT*, 'phyetat0: Le champ <n_snows> est absent' 1035 PRINT *, 'fichier startsisvat non compatible avec sisvatetat0' 1036 ENDIF 1037 1038 CALL get_field("n_ice_top",ispi,found) 1039 CALL get_field("n_ice",iice,found) 1040 CALL get_field("surf_water",rusn,found) 1041 ! IF (.NOT. found) THEN 1042 ! PRINT*, 'phyetat0: Le champ <surf_water> est absent' 1043 ! rusn(:)=0. 1044 ! ENDIF 1045 1046 1047 CALL get_field("to_ice",toic,found) 1048 IF (.NOT. found) THEN 1049 PRINT*, 'phyetat0: Le champ <to_ice> est absent' 1050 toic(:)=0. 1051 ENDIF 1052 1053 1054 1055 DO isn = 1,nsno 1056 IF (isn.LE.99) THEN 1057 WRITE(str2,'(i2.2)') isn 1058 CALL get_field("AGESNOW"//str2, & 1059 agsn(:,isn),found) 1060 ELSE 1061 PRINT*, "Trop de couches" 1062 CALL abort 912 913 !C +--Soil/Snow Layers Indices 914 !C + ^^^^^^^^^^^^^^^^^^^^^^^^ 915 DO isl = -nsol, 0 916 islpSV(isl) = isl + 1 917 islpSV(isl) = min(islpSV(isl), 0) 918 islmSV(isl) = isl - 1 919 islmSV(isl) = max(-nsol, islmSV(isl)) 920 END DO 921 922 DO isn = 1, nsno 923 isnpSV(isn) = isn + 1 924 isnpSV(isn) = min(isnpSV(isn), nsno) 925 END DO 926 927 !C +--Soil Layers Thicknesses 928 !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 929 ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels! 930 !c #kd IF (nsol.gt.4) THEN 931 !c #kd DO isl=-5,-nsol,-1 932 !c #kd dz_dSV(isl)= 1. 933 !c #kd END DO 934 !c #kd END IF 935 ! 936 ! IF (nsol.ne.4) THEN 937 ! DO isl= 0,-nsol,-1 938 ! misl_2 = -mod(isl,2) 939 ! nisl_2 = -isl/2 940 ! dz_dSV(isl)=(((1-misl_2) * 0.001 941 ! . + misl_2 * 0.003) * 10**(nisl_2)) * 4. 942 !C +... dz_dSV(0) = Hapex-Sahel Calibration: 4 mm 943 ! 944 !c +SH dz_dSV(isl)=(((1-misl_2) * 0.001 945 !c +SH. + misl_2 * 0.003) * 10**(nisl_2)) * 1. 946 ! 947 !c #05 dz_dSV(isl)=(((1-misl_2) * 0.001 948 !c #05. + misl_2 * 0.008) * 10**(nisl_2)) * 0.5 949 ! END DO 950 ! dz_dSV(0) = 0.001 951 ! dz_dSV(-1) = dz_dSV(-1) - dz_dSV(0) + 0.004 952 ! END IF 953 954 zz_dSV = 0. 955 DO isl = -nsol, 0 956 dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl))) 957 dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl) 958 dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl) 959 dtz_SV(isl) = dt__SV / dz_dSV(isl) 960 dtz_SV2(isl) = 1. / dz_dSV(isl) 961 dz78SV(isl) = 0.875 * dz_dSV(isl) 962 dz34SV(isl) = 0.750 * dz_dSV(isl) 963 dz_8SV(isl) = 0.125 * dz_dSV(isl) 964 dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl)) & 965 & + 0.750 * dz_dSV(isl) & 966 & + 0.125 * dz_dSV(islpSV(isl)) 967 zz_dSV = zz_dSV + dz_dSV(isl) 968 END DO 969 DO ikl = 1, knon !v 970 dzsnSV(ikl, 0) = dz_dSV(0) 971 END DO 972 973 !C +--Conversion to a 50 m Swab Ocean Discretization 974 !C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 975 OcndSV = 0. 976 DO isl = -nsol, 0 977 OcndSV = OcndSV + dz_dSV(isl) 978 END DO 979 OcndSV = 50. / OcndSV 980 981 982 !C +--Secondary Soil Parameters 983 !C + ------------------------------- 984 985 DO ist = 0, nsot 986 rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6 ! Soil Contrib. to (ro c)_s 987 s1__SV(ist) = bCHdSV(ist) & ! Factor of (eta)**(b+2) 988 & * psidSV(ist) * Ks_dSV(ist) & ! in DR97, Eqn.(3.36) 989 & / (etadSV(ist)**(bCHdSV(ist) + 3.)) ! 990 s2__SV(ist) = Ks_dSV(ist) & ! Factor of (eta)**(2b+3) 991 & / (etadSV(ist)**(2. * bCHdSV(ist) + 3.)) ! in DR97, Eqn.(3.35) 992 993 !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity) 994 !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 995 Psimax = -(log(RHsMin)) / 7.2E-5 ! DR97, Eqn 3.15 Inversion 996 etamSV(ist) = etadSV(ist) & 997 & * (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist))) 998 END DO 999 etamSV(12) = 0. 1000 1001 !C +--Piecewise Hydraulic Conductivity Profiles 1002 !C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1003 DO ist = 0, nsot 1004 1005 d__eta = etadSV(ist) / nkhy 1006 eta__1 = 0. 1007 eta__2 = d__eta 1008 DO ikh = 0, nkhy 1009 Khyd_1 = s2__SV(ist) & ! DR97, Eqn.(3.35) 1010 & * (eta__1 **(2. * bCHdSV(ist) + 3.)) ! 1011 Khyd_2 = s2__SV(ist) &! 1012 & * (eta__2 **(2. * bCHdSV(ist) + 3.)) ! 1013 1014 a_Khyd = (Khyd_2 - Khyd_1) / d__eta ! 1015 b_Khyd = Khyd_1 - a_Khyd * eta__1 ! 1016 !c #WR Khyd_x = a_Khyd*eta__1 +b_Khyd ! 1017 !c #WR Khyd_y = a_Khyd*eta__2 +b_Khyd ! 1018 aKdtSV(ist, ikh) = a_Khyd * dt__SV ! 1019 bKdtSV(ist, ikh) = b_Khyd * dt__SV ! 1020 1021 eta__1 = eta__1 + d__eta 1022 eta__2 = eta__2 + d__eta 1023 END DO 1024 END DO 1025 1026 return 1027 1028 END SUBROUTINE SISVAT_ini 1029 1030 1031 !*************************************************************************** 1032 1033 SUBROUTINE sisvatetat0 (fichnom, ikl2i) 1034 1035 USE dimphy 1036 USE mod_grid_phy_lmdz 1037 USE mod_phys_lmdz_para 1038 1039 USE iostart 1040 USE VAR_SV 1041 USE VARdSV 1042 USE VARxSV 1043 USE VARtSV 1044 USE indice_sol_mod 1045 1046 IMPLICIT none 1047 !====================================================================== 1048 ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009 1049 ! Objet: Lecture du fichier de conditions initiales pour SISVAT 1050 !====================================================================== 1051 include "netcdf.inc" 1052 ! include "indicesol.h" 1053 1054 ! include "dimsoil.h" 1055 include "clesphys.h" 1056 include "thermcell.h" 1057 include "compbl.h" 1058 1059 !====================================================================== 1060 CHARACTER(LEN = *) :: fichnom 1061 1062 INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i 1063 REAL, DIMENSION(klon) :: rlon 1064 REAL, DIMENSION(klon) :: rlat 1065 1066 ! les variables globales ecrites dans le fichier restart 1067 REAL, DIMENSION(klon) :: isno 1068 REAL, DIMENSION(klon) :: ispi 1069 REAL, DIMENSION(klon) :: iice 1070 REAL, DIMENSION(klon) :: rusn 1071 REAL, DIMENSION(klon, nsno) :: isto 1072 1073 REAL, DIMENSION(klon, nsismx) :: Tsis 1074 REAL, DIMENSION(klon, nsismx) :: eta 1075 REAL, DIMENSION(klon, nsismx) :: ro 1076 1077 REAL, DIMENSION(klon, nsno) :: dzsn 1078 REAL, DIMENSION(klon, nsno) :: G1sn 1079 REAL, DIMENSION(klon, nsno) :: G2sn 1080 REAL, DIMENSION(klon, nsno) :: agsn 1081 1082 REAL, DIMENSION(klon) :: toic 1083 1084 INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts 1085 CHARACTER (len = 2) :: str2 1086 LOGICAL :: found 1087 1088 errT = 0 1089 errro = 0 1090 erreta = 0 1091 errdz = 0 1092 snopts = 0 1093 ! Ouvrir le fichier contenant l'etat initial: 1094 1095 CALL open_startphy(fichnom) 1096 1097 ! Lecture des latitudes, longitudes (coordonnees): 1098 1099 CALL get_field("latitude", rlat, found) 1100 CALL get_field("longitude", rlon, found) 1101 1102 CALL get_field("n_snows", isno, found) 1103 IF (.NOT. found) THEN 1104 PRINT*, 'phyetat0: Le champ <n_snows> est absent' 1105 PRINT *, 'fichier startsisvat non compatible avec sisvatetat0' 1063 1106 ENDIF 1064 ENDDO 1065 DO isn = 1,nsno 1066 IF (isn.LE.99) THEN 1067 WRITE(str2,'(i2.2)') isn 1068 CALL get_field("DZSNOW"//str2, & 1069 dzsn(:,isn),found) 1070 ELSE 1071 PRINT*, "Trop de couches" 1072 CALL abort 1107 1108 CALL get_field("n_ice_top", ispi, found) 1109 CALL get_field("n_ice", iice, found) 1110 CALL get_field("surf_water", rusn, found) 1111 1112 1113 CALL get_field("to_ice", toic, found) 1114 IF (.NOT. found) THEN 1115 PRINT*, 'phyetat0: Le champ <to_ice> est absent' 1116 toic(:) = 0. 1073 1117 ENDIF 1074 ENDDO 1075 DO isn = 1,nsno 1076 IF (isn.LE.99) THEN 1077 WRITE(str2,'(i2.2)') isn 1078 CALL get_field("G2SNOW"//str2, & 1079 G2sn(:,isn),found) 1080 ELSE 1081 PRINT*, "Trop de couches" 1082 CALL abort 1083 ENDIF 1084 ENDDO 1085 DO isn = 1,nsno 1086 IF (isn.LE.99) THEN 1087 WRITE(str2,'(i2.2)') isn 1088 CALL get_field("G1SNOW"//str2, & 1089 G1sn(:,isn),found) 1090 ELSE 1091 PRINT*, "Trop de couches" 1092 CALL abort 1093 ENDIF 1094 ENDDO 1095 DO isn = 1,nsismx 1096 IF (isn.LE.99) THEN 1097 WRITE(str2,'(i2.2)') isn 1098 CALL get_field("ETA"//str2, & 1099 eta(:,isn),found) 1100 ELSE 1101 PRINT*, "Trop de couches" 1102 CALL abort 1103 ENDIF 1104 ENDDO 1105 DO isn = 1,nsismx 1106 IF (isn.LE.99) THEN 1107 WRITE(str2,'(i2.2)') isn 1108 CALL get_field("RO"//str2, & 1109 ro(:,isn),found) 1110 ELSE 1111 PRINT*, "Trop de couches" 1112 CALL abort 1113 ENDIF 1114 ENDDO 1115 DO isn = 1,nsismx 1116 IF (isn.LE.99) THEN 1117 WRITE(str2,'(i2.2)') isn 1118 CALL get_field("TSS"//str2, & 1119 Tsis(:,isn),found) 1120 ELSE 1121 PRINT*, "Trop de couches" 1122 CALL abort 1123 ENDIF 1124 ENDDO 1125 DO isn = 1,nsno 1126 IF (isn.LE.99) THEN 1127 WRITE(str2,'(i2.2)') isn 1128 CALL get_field("HISTORY"//str2, & 1129 isto(:,isn),found) 1130 ELSE 1131 PRINT*, "Trop de couches" 1132 CALL abort 1133 ENDIF 1134 ENDDO 1135 write(*,*)'Read ',fichnom,' finished!!' 1136 1137 !********************************************************************************* 1138 ! Compress restart file variables for SISVAT 1139 1140 1141 DO ikl = 1,klon 1142 i = ikl2i(ikl) 1143 IF (i > 0) THEN 1144 isnoSV(ikl) = INT(isno(i)) ! Nb Snow/Ice Lay. 1145 ispiSV(ikl) = INT(ispi(i)) ! Nb Supr.Ice Lay. 1146 iiceSV(ikl) = INT(iice(i)) ! Nb Ice Lay. 1147 1148 1149 DO isl = -nsol,0 1150 ro__SV(ikl,isl) = ro(i,nsno+1-isl) ! 1151 eta_SV(ikl,isl) = eta(i,nsno+1-isl) ! Soil Humidity 1152 !hjp 15/10/2010 1153 IF (eta_SV(ikl,isl) <= 1.e-6) THEN !hj check 1154 eta_SV(ikl,isl) = 1.e-6 1155 ENDIF 1156 TsisSV(ikl,isl) = Tsis(i,nsno+1-isl) ! Soil Temperature 1157 IF (TsisSV(ikl,isl) <= 1.) THEN !hj check 1158 ! errT=errT+1 1159 TsisSV(ikl,isl) = 273.15-0.2 ! Etienne: negative temperature since soil is ice 1160 ENDIF 1161 1162 END DO 1163 write(*,*)'Copy histo', ikl 1164 1165 1166 DO isn = 1,isnoSV(ikl) !nsno 1167 snopts=snopts+1 1168 IF (isto(i,isn) > 10.) THEN !hj check 1169 write(*,*)'Irregular isto',ikl,i,isn,isto(i,isn) 1170 isto(i,isn) = 1. 1171 ENDIF 1172 1173 istoSV(ikl,isn) = INT(isto(i,isn)) ! Snow History 1174 ro__SV(ikl,isn) = ro(i,isn) ! [kg/m3] 1175 eta_SV(ikl,isn) = eta(i,isn) ! [m3/m3] 1176 TsisSV(ikl,isn) = Tsis(i,isn) ! [K] 1177 1178 IF (TsisSV(ikl,isn) <= 1.) THEN !hj check 1179 errT=errT+1 1180 TsisSV(ikl,isn) = TsisSV(ikl,0) 1181 ENDIF 1182 IF (TsisSV(ikl,isn) <= 1.) THEN !hj check 1183 TsisSV(ikl,isn) = 263.15 1184 ENDIF 1185 IF (eta_SV(ikl,isn) < 1.e-9) THEN !hj check 1186 eta_SV(ikl,isn) = 1.e-6 1187 erreta=erreta+1 1188 ENDIF 1189 IF (ro__SV(ikl,isn) <= 10.) THEN !hj check 1190 ro__SV(ikl,isn) = 11. 1191 errro=errro+1 1192 ENDIF 1193 write(*,*)ikl,i,isn,Tsis(i,isn),G1sn(i,isn) 1194 G1snSV(ikl,isn) = G1sn(i,isn) ! [-] [-] 1195 G2snSV(ikl,isn) = G2sn(i,isn) ! [-] [0.0001 m] 1196 dzsnSV(ikl,isn) = dzsn(i,isn) ! [m] 1197 agsnSV(ikl,isn) = agsn(i,isn) ! [day] 1198 END DO 1199 rusnSV(ikl) = rusn(i) ! Surficial Water 1200 toicSV(ikl) = toic(i) ! bilan snow to ice 1201 END IF 1202 END DO 1118 1119 DO isn = 1, nsno 1120 IF (isn.LE.99) THEN 1121 WRITE(str2, '(i2.2)') isn 1122 CALL get_field("AGESNOW" // str2, & 1123 agsn(:, isn), found) 1124 ELSE 1125 PRINT*, "Trop de couches" 1126 CALL abort 1127 ENDIF 1128 ENDDO 1129 DO isn = 1, nsno 1130 IF (isn.LE.99) THEN 1131 WRITE(str2, '(i2.2)') isn 1132 CALL get_field("DZSNOW" // str2, & 1133 dzsn(:, isn), found) 1134 ELSE 1135 PRINT*, "Trop de couches" 1136 CALL abort 1137 ENDIF 1138 ENDDO 1139 DO isn = 1, nsno 1140 IF (isn.LE.99) THEN 1141 WRITE(str2, '(i2.2)') isn 1142 CALL get_field("G2SNOW" // str2, & 1143 G2sn(:, isn), found) 1144 ELSE 1145 PRINT*, "Trop de couches" 1146 CALL abort 1147 ENDIF 1148 ENDDO 1149 DO isn = 1, nsno 1150 IF (isn.LE.99) THEN 1151 WRITE(str2, '(i2.2)') isn 1152 CALL get_field("G1SNOW" // str2, & 1153 G1sn(:, isn), found) 1154 ELSE 1155 PRINT*, "Trop de couches" 1156 CALL abort 1157 ENDIF 1158 ENDDO 1159 DO isn = 1, nsismx 1160 IF (isn.LE.99) THEN 1161 WRITE(str2, '(i2.2)') isn 1162 CALL get_field("ETA" // str2, & 1163 eta(:, isn), found) 1164 ELSE 1165 PRINT*, "Trop de couches" 1166 CALL abort 1167 ENDIF 1168 ENDDO 1169 DO isn = 1, nsismx 1170 IF (isn.LE.99) THEN 1171 WRITE(str2, '(i2.2)') isn 1172 CALL get_field("RO" // str2, & 1173 ro(:, isn), found) 1174 ELSE 1175 PRINT*, "Trop de couches" 1176 CALL abort 1177 ENDIF 1178 ENDDO 1179 DO isn = 1, nsismx 1180 IF (isn.LE.99) THEN 1181 WRITE(str2, '(i2.2)') isn 1182 CALL get_field("TSS" // str2, & 1183 Tsis(:, isn), found) 1184 ELSE 1185 PRINT*, "Trop de couches" 1186 CALL abort 1187 ENDIF 1188 ENDDO 1189 DO isn = 1, nsno 1190 IF (isn.LE.99) THEN 1191 WRITE(str2, '(i2.2)') isn 1192 CALL get_field("HISTORY" // str2, & 1193 isto(:, isn), found) 1194 ELSE 1195 PRINT*, "Trop de couches" 1196 CALL abort 1197 ENDIF 1198 ENDDO 1199 write(*, *)'Read ', fichnom, ' finished!!' 1200 1201 !********************************************************************************* 1202 ! Compress restart file variables for SISVAT 1203 1204 DO ikl = 1, klon 1205 i = ikl2i(ikl) 1206 IF (i > 0) THEN 1207 isnoSV(ikl) = INT(isno(i)) ! Nb Snow/Ice Lay. 1208 ispiSV(ikl) = INT(ispi(i)) ! Nb Supr.Ice Lay. 1209 iiceSV(ikl) = INT(iice(i)) ! Nb Ice Lay. 1210 1211 DO isl = -nsol, 0 1212 ro__SV(ikl, isl) = ro(i, nsno + 1 - isl) ! 1213 eta_SV(ikl, isl) = eta(i, nsno + 1 - isl) ! Soil Humidity 1214 !hjp 15/10/2010 1215 IF (eta_SV(ikl, isl) <= 1.e-6) THEN !hj check 1216 eta_SV(ikl, isl) = 1.e-6 1217 ENDIF 1218 TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl) ! Soil Temperature 1219 IF (TsisSV(ikl, isl) <= 1.) THEN !hj check 1220 ! errT=errT+1 1221 TsisSV(ikl, isl) = 273.15 - 0.2 ! Etienne: negative temperature since soil is ice 1222 ENDIF 1223 1224 END DO 1225 write(*, *)'Copy histo', ikl 1226 1227 DO isn = 1, isnoSV(ikl) !nsno 1228 snopts = snopts + 1 1229 IF (isto(i, isn) > 10.) THEN !hj check 1230 write(*, *)'Irregular isto', ikl, i, isn, isto(i, isn) 1231 isto(i, isn) = 1. 1232 ENDIF 1233 1234 istoSV(ikl, isn) = INT(isto(i, isn)) ! Snow History 1235 ro__SV(ikl, isn) = ro(i, isn) ! [kg/m3] 1236 eta_SV(ikl, isn) = eta(i, isn) ! [m3/m3] 1237 TsisSV(ikl, isn) = Tsis(i, isn) ! [K] 1238 1239 IF (TsisSV(ikl, isn) <= 1.) THEN !hj check 1240 errT = errT + 1 1241 TsisSV(ikl, isn) = TsisSV(ikl, 0) 1242 ENDIF 1243 IF (TsisSV(ikl, isn) <= 1.) THEN !hj check 1244 TsisSV(ikl, isn) = 263.15 1245 ENDIF 1246 IF (eta_SV(ikl, isn) < 1.e-9) THEN !hj check 1247 eta_SV(ikl, isn) = 1.e-6 1248 erreta = erreta + 1 1249 ENDIF 1250 IF (ro__SV(ikl, isn) <= 10.) THEN !hj check 1251 ro__SV(ikl, isn) = 11. 1252 errro = errro + 1 1253 ENDIF 1254 write(*, *)ikl, i, isn, Tsis(i, isn), G1sn(i, isn) 1255 G1snSV(ikl, isn) = G1sn(i, isn) ! [-] [-] 1256 G2snSV(ikl, isn) = G2sn(i, isn) ! [-] [0.0001 m] 1257 dzsnSV(ikl, isn) = dzsn(i, isn) ! [m] 1258 agsnSV(ikl, isn) = agsn(i, isn) ! [day] 1259 END DO 1260 rusnSV(ikl) = rusn(i) ! Surficial Water 1261 toicSV(ikl) = toic(i) ! bilan snow to ice 1262 END IF 1263 END DO 1203 1264 1204 1265 END SUBROUTINE sisvatetat0 1205 1266 1206 1267 1207 1208 1209 !====================================================================== 1210 SUBROUTINE sisvatredem (fichnom,ikl2i,rlon,rlat) 1211 1212 1213 1214 !====================================================================== 1215 ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009 1216 ! Objet: Ecriture de l'etat de redemarrage pour SISVAT 1217 !====================================================================== 1218 USE mod_grid_phy_lmdz 1219 USE mod_phys_lmdz_para 1220 USE iostart 1221 USE VAR_SV 1222 USE VARxSV 1223 USE VARySV !hj tmp 12 03 2010 1224 USE VARtSV 1225 USE indice_sol_mod 1226 USE dimphy 1227 1228 1229 IMPLICIT none 1230 1231 include "netcdf.inc" 1232 ! include "indicesol.h" 1233 ! include "dimsoil.h" 1234 include "clesphys.h" 1235 include "thermcell.h" 1236 include "compbl.h" 1237 1238 !====================================================================== 1239 1240 CHARACTER(LEN=*) :: fichnom 1241 INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i 1242 REAL, DIMENSION(klon), INTENT(IN) :: rlon 1243 REAL, DIMENSION(klon), INTENT(IN) :: rlat 1244 1245 ! les variables globales ecrites dans le fichier restart 1246 REAL, DIMENSION(klon) :: isno 1247 REAL, DIMENSION(klon) :: ispi 1248 REAL, DIMENSION(klon) :: iice 1249 REAL, DIMENSION(klon, nsnowmx) :: isto 1250 1251 REAL, DIMENSION(klon, nsismx) :: Tsis 1252 REAL, DIMENSION(klon, nsismx) :: eta 1253 REAL, DIMENSION(klon, nsnowmx) :: dzsn 1254 REAL, DIMENSION(klon, nsismx) :: ro 1255 REAL, DIMENSION(klon, nsnowmx) :: G1sn 1256 REAL, DIMENSION(klon, nsnowmx) :: G2sn 1257 REAL, DIMENSION(klon, nsnowmx) :: agsn 1258 REAL, DIMENSION(klon) :: IRs 1259 REAL, DIMENSION(klon) :: LMO 1260 REAL, DIMENSION(klon) :: rusn 1261 REAL, DIMENSION(klon) :: toic 1262 REAL, DIMENSION(klon) :: Bufs 1263 REAL, DIMENSION(klon) :: alb1,alb2,alb3 1264 1265 INTEGER isl, ikl, i, isn, ierr 1266 CHARACTER (len=2) :: str2 1267 INTEGER :: pass 1268 1269 isno(:) = 0 1270 ispi(:) = 0 1271 iice(:) = 0 1272 IRs(:) = 0. 1273 LMO(:) = 0. 1274 eta(:,:) = 0. 1275 Tsis(:,:) = 0. 1276 isto(:,:) = 0 1277 ro(:,:) = 0. 1278 G1sn(:,:) = 0. 1279 G2sn(:,:) = 0. 1280 dzsn(:,:) = 0. 1281 agsn(:,:) = 0. 1282 rusn(:) = 0. 1283 toic(:) = 0. 1284 Bufs(:) = 0. 1285 alb1(:) = 0. 1286 alb2(:) = 0. 1287 alb3(:) = 0. 1288 1289 !*************************************************************************** 1290 ! Uncompress SISVAT output variables for storage 1291 1292 1293 print*, 'je rentre dans restart inlandsis' 1294 DO ikl = 1,klon 1295 i = ikl2i(ikl) 1296 IF (i > 0) THEN 1297 isno(i) = 1.*isnoSV(ikl) ! Nb Snow/Ice Lay. 1298 ispi(i) = 1.*ispiSV(ikl) ! Nb Supr.Ice Lay. 1299 iice(i) = 1.*iiceSV(ikl) ! Nb Ice Lay. 1300 1301 ! IRs(i) = IRs_SV(ikl) 1302 ! LMO(i) = LMO_SV(ikl) 1303 1304 1305 DO isl = -nsol,0 ! 1306 eta(i,nsno+1-isl) = eta_SV(ikl,isl) ! Soil Humidity 1307 Tsis(i,nsno+1-isl) = TsisSV(ikl,isl) ! Soil Temperature 1308 ro(i,nsno+1-isl) = ro__SV(ikl,isl) ! [kg/m3] 1309 END DO 1310 1311 1312 DO isn = 1,nsno 1313 isto(i,isn) = 1.*istoSV(ikl,isn) ! Snow History 1314 ro(i,isn) = ro__SV(ikl,isn) ! [kg/m3] 1315 eta(i,isn) = eta_SV(ikl,isn) ! [m3/m3] 1316 Tsis(i,isn) = TsisSV(ikl,isn) ! [K] 1317 G1sn(i,isn) = G1snSV(ikl,isn) ! [-] [-] 1318 G2sn(i,isn) = G2snSV(ikl,isn) ! [-] [0.0001 m] 1319 dzsn(i,isn) = dzsnSV(ikl,isn) ! [m] 1320 agsn(i,isn) = agsnSV(ikl,isn) ! [day] 1321 END DO 1322 rusn(i) = rusnSV(ikl) ! Surficial Water 1323 toic(i) = toicSV(ikl) ! to ice 1324 alb1(i) = alb1sv(ikl) 1325 alb2(i) = alb2sv(ikl) 1326 alb3(i) = alb3sv(ikl) 1327 ! Bufs(i) = BufsSV(ikl) 1328 END IF 1329 END DO 1330 1331 1332 print*, 'je call open_restart' 1333 1334 CALL open_restartphy(fichnom) 1335 1336 print*, 'je sors open_restart' 1337 1338 1339 DO pass = 1, 2 1340 CALL put_field(pass,"longitude", & 1341 "Longitudes de la grille physique",rlon) 1342 CALL put_field(pass,"latitude","Latitudes de la grille physique",rlat) 1343 1344 CALL put_field(pass,"n_snows", "number of snow/ice layers",isno) 1345 CALL put_field(pass,"n_ice_top", "number of top ice layers",ispi) 1346 CALL put_field(pass,"n_ice", "number of ice layers",iice) 1347 CALL put_field(pass,"IR_soil", "Soil IR flux",IRs) 1348 CALL put_field(pass,"LMO", "Monin-Obukhov Scale",LMO) 1349 CALL put_field(pass,"surf_water", "Surficial water",rusn) 1350 CALL put_field(pass,"snow_buffer", "Snow buffer layer",Bufs) 1351 CALL put_field(pass,"alb_1", "albedo sw",alb1) 1352 CALL put_field(pass,"alb_2", "albedo nIR",alb2) 1353 CALL put_field(pass,"alb_3", "albedo fIR",alb3) 1354 CALL put_field(pass,"to_ice", "Snow passed to ice",toic) 1355 1356 1357 1358 DO isn = 1,nsno 1359 IF (isn.LE.99) THEN 1360 WRITE(str2,'(i2.2)') isn 1361 CALL put_field(pass,"AGESNOW"//str2, & 1362 "Age de la neige layer No."//str2, & 1363 agsn(:,isn)) 1364 ELSE 1365 PRINT*, "Trop de couches" 1366 CALL abort 1367 ENDIF 1268 !====================================================================== 1269 SUBROUTINE sisvatredem (fichnom, ikl2i, rlon, rlat) 1270 1271 1272 1273 !====================================================================== 1274 ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009 1275 ! Objet: Ecriture de l'etat de redemarrage pour SISVAT 1276 !====================================================================== 1277 USE mod_grid_phy_lmdz 1278 USE mod_phys_lmdz_para 1279 USE iostart 1280 USE VAR_SV 1281 USE VARxSV 1282 USE VARySV !hj tmp 12 03 2010 1283 USE VARtSV 1284 USE indice_sol_mod 1285 USE dimphy 1286 1287 IMPLICIT none 1288 1289 include "netcdf.inc" 1290 ! include "indicesol.h" 1291 ! include "dimsoil.h" 1292 include "clesphys.h" 1293 include "thermcell.h" 1294 include "compbl.h" 1295 1296 !====================================================================== 1297 1298 CHARACTER(LEN = *) :: fichnom 1299 INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i 1300 REAL, DIMENSION(klon), INTENT(IN) :: rlon 1301 REAL, DIMENSION(klon), INTENT(IN) :: rlat 1302 1303 ! les variables globales ecrites dans le fichier restart 1304 REAL, DIMENSION(klon) :: isno 1305 REAL, DIMENSION(klon) :: ispi 1306 REAL, DIMENSION(klon) :: iice 1307 REAL, DIMENSION(klon, nsnowmx) :: isto 1308 1309 REAL, DIMENSION(klon, nsismx) :: Tsis 1310 REAL, DIMENSION(klon, nsismx) :: eta 1311 REAL, DIMENSION(klon, nsnowmx) :: dzsn 1312 REAL, DIMENSION(klon, nsismx) :: ro 1313 REAL, DIMENSION(klon, nsnowmx) :: G1sn 1314 REAL, DIMENSION(klon, nsnowmx) :: G2sn 1315 REAL, DIMENSION(klon, nsnowmx) :: agsn 1316 REAL, DIMENSION(klon) :: IRs 1317 REAL, DIMENSION(klon) :: LMO 1318 REAL, DIMENSION(klon) :: rusn 1319 REAL, DIMENSION(klon) :: toic 1320 REAL, DIMENSION(klon) :: Bufs 1321 REAL, DIMENSION(klon) :: alb1, alb2, alb3 1322 1323 INTEGER isl, ikl, i, isn, ierr 1324 CHARACTER (len = 2) :: str2 1325 INTEGER :: pass 1326 1327 isno(:) = 0 1328 ispi(:) = 0 1329 iice(:) = 0 1330 IRs(:) = 0. 1331 LMO(:) = 0. 1332 eta(:, :) = 0. 1333 Tsis(:, :) = 0. 1334 isto(:, :) = 0 1335 ro(:, :) = 0. 1336 G1sn(:, :) = 0. 1337 G2sn(:, :) = 0. 1338 dzsn(:, :) = 0. 1339 agsn(:, :) = 0. 1340 rusn(:) = 0. 1341 toic(:) = 0. 1342 Bufs(:) = 0. 1343 alb1(:) = 0. 1344 alb2(:) = 0. 1345 alb3(:) = 0. 1346 1347 !*************************************************************************** 1348 ! Uncompress SISVAT output variables for storage 1349 1350 DO ikl = 1, klon 1351 i = ikl2i(ikl) 1352 IF (i > 0) THEN 1353 isno(i) = 1. * isnoSV(ikl) ! Nb Snow/Ice Lay. 1354 ispi(i) = 1. * ispiSV(ikl) ! Nb Supr.Ice Lay. 1355 iice(i) = 1. * iiceSV(ikl) ! Nb Ice Lay. 1356 1357 ! IRs(i) = IRs_SV(ikl) 1358 ! LMO(i) = LMO_SV(ikl) 1359 1360 DO isl = -nsol, 0 ! 1361 eta(i, nsno + 1 - isl) = eta_SV(ikl, isl) ! Soil Humidity 1362 Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl) ! Soil Temperature 1363 ro(i, nsno + 1 - isl) = ro__SV(ikl, isl) ! [kg/m3] 1364 END DO 1365 1366 DO isn = 1, nsno 1367 isto(i, isn) = 1. * istoSV(ikl, isn) ! Snow History 1368 ro(i, isn) = ro__SV(ikl, isn) ! [kg/m3] 1369 eta(i, isn) = eta_SV(ikl, isn) ! [m3/m3] 1370 Tsis(i, isn) = TsisSV(ikl, isn) ! [K] 1371 G1sn(i, isn) = G1snSV(ikl, isn) ! [-] [-] 1372 G2sn(i, isn) = G2snSV(ikl, isn) ! [-] [0.0001 m] 1373 dzsn(i, isn) = dzsnSV(ikl, isn) ! [m] 1374 agsn(i, isn) = agsnSV(ikl, isn) ! [day] 1375 END DO 1376 rusn(i) = rusnSV(ikl) ! Surficial Water 1377 toic(i) = toicSV(ikl) ! to ice 1378 alb1(i) = alb1sv(ikl) 1379 alb2(i) = alb2sv(ikl) 1380 alb3(i) = alb3sv(ikl) 1381 ! Bufs(i) = BufsSV(ikl) 1382 END IF 1383 END DO 1384 1385 CALL open_restartphy(fichnom) 1386 1387 DO pass = 1, 2 1388 CALL put_field(pass, "longitude", & 1389 "Longitudes de la grille physique", rlon) 1390 CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat) 1391 1392 CALL put_field(pass, "n_snows", "number of snow/ice layers", isno) 1393 CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi) 1394 CALL put_field(pass, "n_ice", "number of ice layers", iice) 1395 CALL put_field(pass, "IR_soil", "Soil IR flux", IRs) 1396 CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO) 1397 CALL put_field(pass, "surf_water", "Surficial water", rusn) 1398 CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs) 1399 CALL put_field(pass, "alb_1", "albedo sw", alb1) 1400 CALL put_field(pass, "alb_2", "albedo nIR", alb2) 1401 CALL put_field(pass, "alb_3", "albedo fIR", alb3) 1402 CALL put_field(pass, "to_ice", "Snow passed to ice", toic) 1403 1404 DO isn = 1, nsno 1405 IF (isn.LE.99) THEN 1406 WRITE(str2, '(i2.2)') isn 1407 CALL put_field(pass, "AGESNOW" // str2, & 1408 "Age de la neige layer No." // str2, & 1409 agsn(:, isn)) 1410 ELSE 1411 PRINT*, "Trop de couches" 1412 CALL abort 1413 ENDIF 1414 ENDDO 1415 DO isn = 1, nsno 1416 IF (isn.LE.99) THEN 1417 WRITE(str2, '(i2.2)') isn 1418 CALL put_field(pass, "DZSNOW" // str2, & 1419 "Snow/ice thickness layer No." // str2, & 1420 dzsn(:, isn)) 1421 ELSE 1422 PRINT*, "Trop de couches" 1423 CALL abort 1424 ENDIF 1425 ENDDO 1426 DO isn = 1, nsno 1427 IF (isn.LE.99) THEN 1428 WRITE(str2, '(i2.2)') isn 1429 CALL put_field(pass, "G2SNOW" // str2, & 1430 "Snow Property 2, layer No." // str2, & 1431 G2sn(:, isn)) 1432 ELSE 1433 PRINT*, "Trop de couches" 1434 CALL abort 1435 ENDIF 1436 ENDDO 1437 DO isn = 1, nsno 1438 IF (isn.LE.99) THEN 1439 WRITE(str2, '(i2.2)') isn 1440 CALL put_field(pass, "G1SNOW" // str2, & 1441 "Snow Property 1, layer No." // str2, & 1442 G1sn(:, isn)) 1443 ELSE 1444 PRINT*, "Trop de couches" 1445 CALL abort 1446 ENDIF 1447 ENDDO 1448 DO isn = 1, nsismx 1449 IF (isn.LE.99) THEN 1450 WRITE(str2, '(i2.2)') isn 1451 CALL put_field(pass, "ETA" // str2, & 1452 "Soil/snow water content layer No." // str2, & 1453 eta(:, isn)) 1454 ELSE 1455 PRINT*, "Trop de couches" 1456 CALL abort 1457 ENDIF 1458 ENDDO 1459 DO isn = 1, nsismx !nsno 1460 IF (isn.LE.99) THEN 1461 WRITE(str2, '(i2.2)') isn 1462 CALL put_field(pass, "RO" // str2, & 1463 "Snow density layer No." // str2, & 1464 ro(:, isn)) 1465 ELSE 1466 PRINT*, "Trop de couches" 1467 CALL abort 1468 ENDIF 1469 ENDDO 1470 DO isn = 1, nsismx 1471 IF (isn.LE.99) THEN 1472 WRITE(str2, '(i2.2)') isn 1473 CALL put_field(pass, "TSS" // str2, & 1474 "Soil/snow temperature layer No." // str2, & 1475 Tsis(:, isn)) 1476 ELSE 1477 PRINT*, "Trop de couches" 1478 CALL abort 1479 ENDIF 1480 ENDDO 1481 DO isn = 1, nsno 1482 IF (isn.LE.99) THEN 1483 WRITE(str2, '(i2.2)') isn 1484 CALL put_field(pass, "HISTORY" // str2, & 1485 "Snow history layer No." // str2, & 1486 isto(:, isn)) 1487 ELSE 1488 PRINT*, "Trop de couches" 1489 CALL abort 1490 ENDIF 1491 ENDDO 1492 1493 CALL enddef_restartphy 1368 1494 ENDDO 1369 DO isn = 1,nsno 1370 IF (isn.LE.99) THEN 1371 WRITE(str2,'(i2.2)') isn 1372 CALL put_field(pass,"DZSNOW"//str2, & 1373 "Snow/ice thickness layer No."//str2, & 1374 dzsn(:,isn)) 1375 ELSE 1376 PRINT*, "Trop de couches" 1377 CALL abort 1378 ENDIF 1379 ENDDO 1380 DO isn = 1,nsno 1381 IF (isn.LE.99) THEN 1382 WRITE(str2,'(i2.2)') isn 1383 CALL put_field(pass,"G2SNOW"//str2, & 1384 "Snow Property 2, layer No."//str2, & 1385 G2sn(:,isn)) 1386 ELSE 1387 PRINT*, "Trop de couches" 1388 CALL abort 1389 ENDIF 1390 ENDDO 1391 DO isn = 1,nsno 1392 IF (isn.LE.99) THEN 1393 WRITE(str2,'(i2.2)') isn 1394 CALL put_field(pass,"G1SNOW"//str2, & 1395 "Snow Property 1, layer No."//str2, & 1396 G1sn(:,isn)) 1397 ELSE 1398 PRINT*, "Trop de couches" 1399 CALL abort 1400 ENDIF 1401 ENDDO 1402 DO isn = 1,nsismx 1403 IF (isn.LE.99) THEN 1404 WRITE(str2,'(i2.2)') isn 1405 CALL put_field(pass,"ETA"//str2, & 1406 "Soil/snow water content layer No."//str2, & 1407 eta(:,isn)) 1408 ELSE 1409 PRINT*, "Trop de couches" 1410 CALL abort 1411 ENDIF 1412 ENDDO 1413 DO isn = 1,nsismx !nsno 1414 IF (isn.LE.99) THEN 1415 WRITE(str2,'(i2.2)') isn 1416 CALL put_field(pass,"RO"//str2, & 1417 "Snow density layer No."//str2, & 1418 ro(:,isn)) 1419 ELSE 1420 PRINT*, "Trop de couches" 1421 CALL abort 1422 ENDIF 1423 ENDDO 1424 DO isn = 1,nsismx 1425 IF (isn.LE.99) THEN 1426 WRITE(str2,'(i2.2)') isn 1427 CALL put_field(pass,"TSS"//str2, & 1428 "Soil/snow temperature layer No."//str2, & 1429 Tsis(:,isn)) 1430 ELSE 1431 PRINT*, "Trop de couches" 1432 CALL abort 1433 ENDIF 1434 ENDDO 1435 DO isn = 1,nsno 1436 IF (isn.LE.99) THEN 1437 WRITE(str2,'(i2.2)') isn 1438 CALL put_field(pass,"HISTORY"//str2, & 1439 "Snow history layer No."//str2, & 1440 isto(:,isn)) 1441 ELSE 1442 PRINT*, "Trop de couches" 1443 CALL abort 1444 ENDIF 1445 ENDDO 1446 1447 CALL enddef_restartphy 1448 ENDDO 1449 CALL close_restartphy 1450 1451 1452 END SUBROUTINE sisvatredem 1495 CALL close_restartphy 1496 1497 END SUBROUTINE sisvatredem 1453 1498 1454 1499 END MODULE surf_inlandsis_mod
Note: See TracChangeset
for help on using the changeset viewer.
