Changeset 3900
- Timestamp:
- May 17, 2021, 3:35:58 PM (4 years ago)
- Location:
- LMDZ6/trunk/libf/phylmd
- Files:
-
- 20 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
r3876 r3900 27 27 USE phys_cal_mod 28 28 USE carbon_cycle_mod, ONLY: carbon_cycle_tr, carbon_cycle_cpl, carbon_cycle_rad, level_coupling_esm 29 USE carbon_cycle_mod, ONLY: read_fco2_ocean_cor,var_fco2_ocean_cor30 USE carbon_cycle_mod, ONLY: read_fco2_land_cor,var_fco2_land_cor31 29 USE mod_grid_phy_lmdz, ONLY: klon_glo 32 30 USE print_control_mod, ONLY: lunout … … 90 88 CHARACTER (len = 8), SAVE :: aer_type_omp 91 89 INTEGER, SAVE :: landice_opt_omp 92 INTEGER, SAVE :: n_dtis_omp 93 INTEGER, SAVE :: iflag_tsurf_inlandsis_omp 94 INTEGER, SAVE :: iflag_albzenith_omp 95 LOGICAL, SAVE :: SnoMod_omp,BloMod_omp,ok_outfor_omp 90 INTEGER, SAVE :: iflag_tsurf_inlandsis_omp,iflag_temp_inlandsis_omp 91 INTEGER, SAVE :: iflag_albcalc_omp,iflag_z0m_snow_omp 92 LOGICAL, SAVE :: SnoMod_omp,BloMod_omp,ok_outfor_omp,ok_zsn_ii_omp 93 LOGICAL, SAVE :: discret_xf_omp,opt_runoff_ac_omp 94 LOGICAL, SAVE :: is_ok_slush_omp,is_ok_z0h_rn_omp,is_ok_density_kotlyakov_omp 95 REAL, SAVE :: prescribed_z0m_snow_omp,correc_alb_omp 96 REAL, SAVE :: buf_sph_pol_omp,buf_siz_pol_omp 96 97 LOGICAL, SAVE :: ok_newmicro_omp 97 98 LOGICAL, SAVE :: ok_all_xml_omp … … 239 240 LOGICAL, SAVE :: carbon_cycle_rad_omp 240 241 INTEGER, SAVE :: level_coupling_esm_omp 241 LOGICAL, SAVE :: read_fco2_ocean_cor_omp242 REAL, SAVE :: var_fco2_ocean_cor_omp243 LOGICAL, SAVE :: read_fco2_land_cor_omp244 REAL, SAVE :: var_fco2_land_cor_omp245 242 LOGICAL, SAVE :: adjust_tropopause_omp 246 243 LOGICAL, SAVE :: ok_daily_climoz_omp … … 340 337 341 338 !Etienne 339 !Config Key = iflag_temp_inlandsis 340 !Config Desc = which method to calculate temp within the soil in INLANDSIS 341 !Config Def = 0 342 iflag_temp_inlandsis_omp = 0 343 CALL getin('iflag_temp_inlandsis', iflag_temp_inlandsis_omp) 344 345 !Etienne 342 346 !Config Key = iflag_tsurf_inlandsis 343 347 !Config Desc = which method to calculate tsurf in INLANDSIS 344 348 !Config Def = 0 345 iflag_tsurf_inlandsis_omp = 0349 iflag_tsurf_inlandsis_omp = 1 346 350 CALL getin('iflag_tsurf_inlandsis', iflag_tsurf_inlandsis_omp) 347 351 352 348 353 !Etienne 349 !Config Key = iflag_albzenith 350 !Config Desc = method to account for albedo sensitivity to solar zenith angle 351 !Config Def = 0 352 iflag_albzenith_omp = 0 353 CALL getin('iflag_albzenith', iflag_albzenith_omp) 354 355 !Etienne 356 !Config Key = n_dtis 357 !Config Desc = number of subtimesteps for INLANDSIS 358 !Config Def = 1 359 n_dtis_omp = 1 360 CALL getin('n_dtis', n_dtis_omp) 354 !Config Key = iflag_albcalc 355 !Config Desc = method to calculate snow albedo in INLANDSIS 356 !Config Def = 0 357 iflag_albcalc_omp = 0 358 CALL getin('iflag_albcalc', iflag_albcalc_omp) 359 361 360 362 361 !Etienne 363 362 !Config Key = SnoMod 364 363 !Config Desc = activation of snow modules in inlandsis 365 !Config Def = 1364 !Config Def = .TRUE. 366 365 SnoMod_omp = .TRUE. 367 366 CALL getin('SnoMod', SnoMod_omp) … … 370 369 !Config Key = BloMod 371 370 !Config Desc = activation of blowing snow in inlandsis 372 !Config Def = 1371 !Config Def = .FALSE. 373 372 BloMod_omp = .FALSE. 374 373 CALL getin('BloMod', BloMod_omp) … … 377 376 !Config Key = ok_outfor 378 377 !Config Desc = activation of output ascii file in inlandsis 379 !Config Def = 1378 !Config Def = .FALSE. 380 379 ok_outfor_omp = .FALSE. 381 380 CALL getin('ok_outfor', ok_outfor_omp) 381 382 383 !Etienne 384 !Config Key = ok_sn_ii 385 !Config Desc = activation of ice/snow detection 386 !Config Def = .TRUE. 387 ok_zsn_ii = .TRUE. 388 CALL getin('ok_zsn_ii', ok_zsn_ii_omp) 389 390 391 !Etienne 392 !Config Key = discret_xf 393 !Config Desc = snow discretization following XF 394 !Config Def = .TRUE. 395 discret_xf = .TRUE. 396 CALL getin('discret_xf', discret_xf_omp) 397 398 399 !Etienne 400 !Config Key = is_ok_slush 401 !Config Desc = activation of the slush option 402 !Config Def = .TRUE. 403 is_ok_slush_omp = .TRUE. 404 CALL getin('is_ok_slush', is_ok_slush_omp) 405 406 !Etienne 407 !Config Key = opt_runoff_ac 408 !Config Desc = option runoff AC 409 !Config Def = .TRUE. 410 opt_runoff_ac_omp = .TRUE. 411 CALL getin('opt_runoff_ac', opt_runoff_ac_omp) 412 413 !Etienne 414 !Config Key = is_ok_z0h_rn 415 !Config Desc = z0h calculation following RN method 416 !Config Def = .TRUE. 417 is_ok_z0h_rn = .TRUE. 418 CALL getin('is_ok_z0h_rn', is_ok_z0h_rn_omp) 419 420 421 !Etienne 422 !Config Key = is_ok_density_kotlyakov 423 !Config Desc = snow density calculation following kotlyakov 424 !Config Def = .FALSE. 425 is_ok_density_kotlyakov = .FALSE. 426 CALL getin('is_ok_density_kotlyakov', is_ok_density_kotlyakov_omp) 427 428 429 !Etienne 430 !Config Key = prescribed_z0m_snow 431 !Config Desc = prescribed snow z0m 432 !Config Def = 0.005 433 prescribed_z0m_snow = 0.005 434 CALL getin('prescribed_z0m_snow', prescribed_z0m_snow_omp) 435 436 437 !Etienne 438 !Config Key = iflag_z0m_snow 439 !Config Desc = method to calculate snow z0m 440 !Config Def = 0 441 iflag_z0m_snow = 0 442 CALL getin('iflag_z0m_snow', iflag_z0m_snow_omp) 443 444 445 !Etienne 446 !Config Key = correc_alb 447 !Config Desc = correction term for albedo 448 !Config Def = 1.01 449 correc_alb_omp=1.01 450 CALL getin('correc_alb', correc_alb_omp) 451 452 453 !Etienne 454 !Config Key = buf_sph_pol 455 !Config Desc = sphericity of buffer layer in polar regions 456 !Config Def = 99. 457 buf_sph_pol_omp=99. 458 CALL getin('buf_sph_pol', buf_sph_pol_omp) 459 460 !Etienne 461 !Config Key = buf_siz_pol 462 !Config Desc = grain size of buffer layer in polar regions in e-4m 463 !Config Def = 4. 464 buf_siz_pol_omp=4. 465 CALL getin('buf_siz_pol', buf_siz_pol_omp) 382 466 383 467 !================================================================== … … 1140 1224 CALL getin('coefw_cld_cv',coefw_cld_cv_omp) 1141 1225 1226 1227 1228 1142 1229 ! 1143 1230 !Config Key = iflag_pdf … … 1148 1235 iflag_pdf_omp = 0 1149 1236 CALL getin('iflag_pdf',iflag_pdf_omp) 1150 1151 1237 ! 1152 1238 !Config Key = fact_cldcon … … 1175 1261 ok_newmicro_omp = .TRUE. 1176 1262 CALL getin('ok_newmicro',ok_newmicro_omp) 1177 1178 1263 ! 1179 1264 !Config Key = ratqsbas … … 1184 1269 ratqsbas_omp = 0.01 1185 1270 CALL getin('ratqsbas',ratqsbas_omp) 1186 1187 1271 ! 1188 1272 !Config Key = ratqshaut … … 2218 2302 CALL getin('carbon_cycle_rad',carbon_cycle_rad_omp) 2219 2303 2220 read_fco2_ocean_cor_omp=.FALSE. 2221 CALL getin('read_fco2_ocean_cor',read_fco2_ocean_cor_omp) 2222 2223 var_fco2_ocean_cor_omp=0. ! default value 2224 CALL getin('var_fco2_ocean_cor',var_fco2_ocean_cor_omp) 2225 2226 read_fco2_land_cor_omp=.FALSE. 2227 CALL getin('read_fco2_land_cor',read_fco2_land_cor_omp) 2228 2229 var_fco2_land_cor_omp=0. ! default value 2230 CALL getin('var_fco2_land_cor',var_fco2_land_cor_omp) 2231 2304 ! >> PC 2232 2305 ! level_coupling_esm : level of coupling of the biogeochemical fields between LMDZ, ORCHIDEE and NEMO 2233 2306 ! Definitions of level_coupling_esm in physiq.def … … 2242 2315 level_coupling_esm_omp=0 ! default value 2243 2316 CALL getin('level_coupling_esm',level_coupling_esm_omp) 2317 ! << PC 2244 2318 2245 2319 !$OMP END MASTER … … 2358 2432 ok_veget=.FALSE. 2359 2433 ENDIF 2360 ! SISVAT andINLANDSIS2434 ! INLANDSIS 2361 2435 !================================================= 2362 2436 landice_opt = landice_opt_omp 2363 2437 iflag_tsurf_inlandsis = iflag_tsurf_inlandsis_omp 2364 iflag_ albzenith = iflag_albzenith_omp2365 n_dtis=n_dtis_omp2438 iflag_temp_inlandsis = iflag_temp_inlandsis_omp 2439 iflag_albcalc = iflag_albcalc_omp 2366 2440 SnoMod=SnoMod_omp 2367 2441 BloMod=BloMod_omp 2368 2442 ok_outfor=ok_outfor_omp 2443 is_ok_slush=is_ok_slush_omp 2444 opt_runoff_ac=opt_runoff_ac_omp 2445 is_ok_z0h_rn=is_ok_z0h_rn_omp 2446 is_ok_density_kotlyakov=is_ok_density_kotlyakov_omp 2447 prescribed_z0m_snow=prescribed_z0m_snow_omp 2448 correc_alb=correc_alb_omp 2449 iflag_z0m_snow=iflag_z0m_snow_omp 2450 ok_zsn_ii=ok_zsn_ii_omp 2451 discret_xf=discret_xf_omp 2452 buf_sph_pol=buf_sph_pol_omp 2453 buf_siz_pol=buf_siz_pol_omp 2369 2454 !================================================= 2370 2455 ok_all_xml = ok_all_xml_omp … … 2507 2592 carbon_cycle_rad = carbon_cycle_rad_omp 2508 2593 level_coupling_esm = level_coupling_esm_omp 2509 read_fco2_ocean_cor = read_fco2_ocean_cor_omp2510 var_fco2_ocean_cor = var_fco2_ocean_cor_omp2511 read_fco2_land_cor = read_fco2_land_cor_omp2512 var_fco2_land_cor = var_fco2_land_cor_omp2513 2594 2514 2595 ! Test of coherence between type_ocean and version_ocean … … 2831 2912 WRITE(lunout,*) ' level_coupling_esm = ', level_coupling_esm 2832 2913 WRITE(lunout,*) ' iflag_tsurf_inlandsis = ', iflag_tsurf_inlandsis 2833 WRITE(lunout,*) ' iflag_ albzenith = ', iflag_albzenith2834 WRITE(lunout,*) ' n_dtis = ', n_dtis2914 WRITE(lunout,*) ' iflag_temp_inlandsis = ', iflag_temp_inlandsis 2915 WRITE(lunout,*) ' iflag_albcalc = ', iflag_albcalc 2835 2916 WRITE(lunout,*) ' SnoMod = ', SnoMod 2836 2917 WRITE(lunout,*) ' BloMod = ', BloMod 2837 2918 WRITE(lunout,*) ' ok_outfor = ', ok_outfor 2838 2919 WRITE(lunout,*) ' is_ok_slush = ', is_ok_slush 2920 WRITE(lunout,*) ' opt_runoff_ac = ', opt_runoff_ac 2921 WRITE(lunout,*) ' is_ok_z0h_rn = ', is_ok_z0h_rn 2922 WRITE(lunout,*) ' is_ok_density_kotlyakov = ', is_ok_density_kotlyakov 2923 WRITE(lunout,*) ' prescribed_z0m_snow = ', prescribed_z0m_snow 2924 WRITE(lunout,*) ' iflag_z0m_snow = ', iflag_z0m_snow 2925 WRITE(lunout,*) ' ok_zsn_ii = ', ok_zsn_ii 2926 WRITE(lunout,*) ' discret_xf = ', discret_xf 2927 WRITE(lunout,*) ' correc_alb= ', correc_alb 2928 WRITE(lunout,*) ' buf_sph_pol = ', buf_sph_pol 2929 WRITE(lunout,*) ' buf_siz_pol= ', buf_siz_pol 2839 2930 2840 2931 !$OMP END MASTER -
LMDZ6/trunk/libf/phylmd/create_etat0_unstruct.F90
r3535 r3900 209 209 z0m(:,is_oce) = rugmer(:) 210 210 211 z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)212 z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)211 z0m(:,is_ter) = 0.01 ! MAX(1.0e-05,zstd(:)*zsig(:)/2.0) 212 z0m(:,is_lic) = 0.001 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0) 213 213 214 214 z0m(:,is_sic) = 0.001 -
LMDZ6/trunk/libf/phylmd/dimsoil.h
r3792 r3900 8 8 9 9 INTEGER nsnowmx 10 PARAMETER (nsnowmx=3 5)10 PARAMETER (nsnowmx=30) 11 11 12 12 INTEGER nsismx 13 PARAMETER (nsismx=4 6)13 PARAMETER (nsismx=41) 14 14 15 15 ! nsismx should be equal to nsoilmx+nsnowmx -
LMDZ6/trunk/libf/phylmd/fonte_neige_mod.F90
r3102 r3900 28 28 REAL, PRIVATE :: tau_calv 29 29 !$OMP THREADPRIVATE(tau_calv) 30 REAL, ALLOCATABLE, DIMENSION(:,:) , PRIVATE:: ffonte_global30 REAL, ALLOCATABLE, DIMENSION(:,:) :: ffonte_global 31 31 !$OMP THREADPRIVATE(ffonte_global) 32 REAL, ALLOCATABLE, DIMENSION(:,:) , PRIVATE:: fqfonte_global32 REAL, ALLOCATABLE, DIMENSION(:,:) :: fqfonte_global 33 33 !$OMP THREADPRIVATE(fqfonte_global) 34 REAL, ALLOCATABLE, DIMENSION(:,:) , PRIVATE:: fqcalving_global34 REAL, ALLOCATABLE, DIMENSION(:,:) :: fqcalving_global 35 35 !$OMP THREADPRIVATE(fqcalving_global) 36 REAL, ALLOCATABLE, DIMENSION(:) , PRIVATE:: runofflic_global36 REAL, ALLOCATABLE, DIMENSION(:) :: runofflic_global 37 37 !$OMP THREADPRIVATE(runofflic_global) 38 38 -
LMDZ6/trunk/libf/phylmd/inlandsis/VARphy.F90
r3792 r3900 26 26 INTEGER, PARAMETER :: iun=1 27 27 REAL, PARAMETER :: zer0 = 0.0e+0, half = 0.5e+0, un_1 = 1.0e+0, & 28 & eps6 = 1.0e-6, R_1000=1.e3 28 & eps6 = 1.0e-6, R_1000=1.e3 29 29 REAL, PARAMETER :: zero = 0.0e+0, demi = 0.5e+0, unun = 1.0e+0, & 30 30 & epsi = 1.0e-6, eps9 = 1.0e-9 … … 91 91 ! A1.6 Turbulent and molecular diffusion 92 92 !---------------------------------------- 93 REAL, PARAMETER :: A_MolV = 1.35e-5, vonKrm = 0.40e0 93 REAL, PARAMETER :: A_MolV = 1.35e-5, vonKrm = 0.40e0, r_turb=3.0 94 REAL, PARAMETER :: A_turb=5.8, akmol=1.35e-5 94 95 !C +... A_MolV: Air Viscosity = 1.35d-5 m2/s 95 96 !C + vonKrm: von Karman constant = 0.4 96 97 !C + r_turb: Turbulent Diffusivities Ratio K*/Km 98 !C + A_turb: Stability Coefficient Moment 99 !C + Air Viscosity = 1.35d-5 m2/s 100 101 97 102 98 103 END MODULE VARphy -
LMDZ6/trunk/libf/phylmd/inlandsis/VARtSV.F90
r3792 r3900 41 41 42 42 SUBROUTINE INIT_VARtSV 43 43 44 IMPLICIT NONE 44 45 46 INTEGER ikl 47 48 49 50 51 52 53 45 54 ALLOCATE(toicSV(klonv)) 46 55 … … 59 68 ALLOCATE(rsolSV(klonv)) ! Radiation balance surface 60 69 70 DO ikl=1,klonv 71 72 toicSV(ikl) = 0. 73 dz1_SV(ikl,:) = 0. 74 dz2_SV(ikl,:) = 0. 75 Tsf_SV(ikl) = 0. 76 TsfnSV(ikl) = 0. 77 AcoHSV(ikl) = 0. 78 BcoHSV(ikl) = 0. 79 AcoQSV(ikl) = 0. 80 ps__SV(ikl) = 0. 81 p1l_SV(ikl) = 0. 82 cdH_SV(ikl) = 0. 83 cdM_SV(ikl) = 0. 84 rsolSV(ikl) = 0. 85 END DO 86 87 88 61 89 END SUBROUTINE INIT_VARtSV 62 90 -
LMDZ6/trunk/libf/phylmd/inlandsis/VARxSV.F90
r3792 r3900 67 67 REAL, DIMENSION(:),ALLOCATABLE,SAVE :: QaT_SV ! SBL Top Specific Humidity 68 68 !$OMP THREADPRIVATE(QaT_SV) 69 REAL, DIMENSION(:),ALLOCATABLE,SAVE :: QsT_SV ! SBL Top Specific Humidity 70 !$OMP THREADPRIVATE(QsT_SV) 69 71 REAL, DIMENSION(:),ALLOCATABLE,SAVE :: dQa_SV ! SBL Flux Limitation of Qa 70 72 !$OMP THREADPRIVATE(dQa_SV) … … 78 80 79 81 80 REAL,SAVE :: zSBLSV ! SBL Height (Initial Value)81 !$OMP THREADPRIVATE(zSBLSV)82 82 REAL,SAVE :: dt__SV ! Time Step 83 83 !$OMP THREADPRIVATE(dt__SV) … … 160 160 REAL,ALLOCATABLE,SAVE :: agsnSV(:,:) ! Snow Age 161 161 !$OMP THREADPRIVATE(agsnSV) 162 REAL,ALLOCATABLE,SAVE :: DOPsnSV(:,:) ! Snow optical diameter [m] 163 !$OMP THREADPRIVATE(DOPsnSV) 162 164 REAL, DIMENSION(:),ALLOCATABLE,SAVE :: BufsSV ! Snow Buffer Layer 163 165 !$OMP THREADPRIVATE(BufsSV) … … 260 262 ALLOCATE(dLdTSV(klonv)) ! Latent Heat Flux T Derivat. 261 263 ALLOCATE(rhT_SV(klonv)) ! SBL Top Air Density 262 ALLOCATE(QaT_SV(klonv)) ! SBL Top Specific Humidity 264 ALLOCATE(QaT_SV(klonv)) ! SBL Top Specific Humidity 265 ALLOCATE(QsT_SV(klonv)) ! surface Specific Humidity 263 266 ALLOCATE(dQa_SV(klonv)) ! SBL Flux Limitation of Qa 264 267 ALLOCATE(qsnoSV(klonv)) ! SBL Mean Snow Content … … 309 312 ALLOCATE(dzsnSV(klonv, 0:nsno)) ! Snow Layer Thickness 310 313 ALLOCATE(agsnSV(klonv, 0:nsno)) ! Snow Age 314 ALLOCATE(DOPsnSV(klonv, 0:nsno)) ! Snow Optical diameter 311 315 ALLOCATE(BufsSV(klonv)) ! Snow Buffer Layer 312 316 ALLOCATE(rusnSV(klonv)) ! Surficial Water … … 339 343 340 344 DO ikl=1,klonv 341 LSmask(ikl) = 0 342 isotSV(ikl) = 0 343 iWaFSV(ikl) = 0 344 isnoSV(ikl) = 0 345 ispiSV(ikl) = 0 346 iiceSV(ikl) = 0 347 istoSV(ikl,:) = 0 348 ii__SV(ikl) = 0 349 jj__SV(ikl) = 0 350 nn__SV(ikl) = 0 345 346 347 isnoSV(ikl) =0. 348 ispiSV(ikl) =0. 349 iiceSV(ikl) =0. 350 istoSV(ikl,:)=0. 351 alb_SV(ikl) =0. 352 emi_SV(ikl) =0. 353 IRs_SV(ikl) =0. 354 LMO_SV(ikl) =0. 355 us__SV(ikl) =0. 356 uts_SV(ikl) =0. 357 cutsSV(ikl) =0. 358 uqs_SV(ikl) =0. 359 uss_SV(ikl) =0. 360 usthSV(ikl) =0. 361 rCDmSV(ikl) =0. 362 rCDhSV(ikl) =0. 363 Z0m_SV(ikl) =0. 364 Z0mmSV(ikl) =0. 365 Z0mnSV(ikl) =0. 366 Z0roSV(ikl) =0. 367 Z0SaSV(ikl) =0. 368 Z0e_SV(ikl) =0. 369 Z0emSV(ikl) =0. 370 Z0enSV(ikl) =0. 371 Z0h_SV(ikl) =0. 372 Z0hmSV(ikl) =0. 373 Z0hnSV(ikl) =0. 374 375 376 TsisSV(ikl,:) =0. 377 ro__SV(ikl,:) =0. 378 eta_SV(ikl,:) =0. 379 G1snSV(ikl,:) =0. 380 G2snSV(ikl,:) =0. 381 dzsnSV(ikl,:) =0. 382 agsnSV(ikl,:) =0. 383 DOPsnSV(ikl,:) =0. 384 BufsSV(ikl) =0. 385 rusnSV(ikl) =0. 386 SWf_SV(ikl) =0. 387 SWS_SV(ikl) =0. 388 HFraSV(ikl) =0. 389 390 zWE_SV(ikl) =0. 391 zWEcSV(ikl) =0. 392 wem_SV(ikl) =0. 393 wer_SV(ikl) =0. 394 wes_SV(ikl) =0. 395 zn4_SV(ikl) =0. 396 zn5_SV(ikl) =0. 397 398 399 ii__SV(ikl) =0. 400 jj__SV(ikl) =0. 401 nn__SV(ikl) =0. 402 403 IRu_SV(ikl) =0. 404 hSalSV(ikl) =0. 405 qSalSV(ikl) =0. 406 RnofSV(ikl) =0. 407 RuofSV(ikl,:) =0. 408 409 410 411 351 412 END DO 352 413 END SUBROUTINE INIT_VARxSV -
LMDZ6/trunk/libf/phylmd/inlandsis/VARySV.F90
r3792 r3900 22 22 REAL, DIMENSION(:),SAVE,ALLOCATABLE :: alb3sv ! Surface Albedo FIR 23 23 !$OMP THREADPRIVATE(alb3sv) 24 24 REAL, DIMENSION(:,:),SAVE,ALLOCATABLE :: alb6sv ! 6 band-albedo 25 !$OMP THREADPRIVATE(alb6sv) 25 26 REAL, DIMENSION(:),SAVE,ALLOCATABLE :: albssv ! Soil Albedo [-] 26 27 !$OMP THREADPRIVATE(albssv) … … 83 84 ALLOCATE(alb2sv(klonv)) ! Surface Albedo NIR 84 85 ALLOCATE(alb3sv(klonv)) ! Surface Albedo FIR 86 ALLOCATE(alb6sv(klonv,6))! 6-band Albedo 85 87 86 88 ! … … 110 112 111 113 DO ikl=1,klonv 112 NLaysv(ikl) = 0 113 i_thin(ikl) = 0 114 LIndsv(ikl) = 0 114 115 NLaysv(ikl) =0. 116 i_thin(ikl) =0. 117 LIndsv(ikl) =0. 118 albisv(ikl) =0. 119 alb1sv(ikl) =0. 120 alb2sv(ikl) =0. 121 alb3sv(ikl) =0. 122 alb6sv(ikl,:)=0. 123 albssv(ikl) =0. 124 SoSosv(ikl) =0. 125 Eso_sv(ikl) =0. 126 HSv_sv(ikl) =0. 127 HLv_sv(ikl) =0. 128 HSs_sv(ikl) =0. 129 HLs_sv(ikl) =0. 130 sqrCm0(ikl) =0. 131 sqrCh0(ikl) =0. 132 Lx_H2O(ikl) =0. 133 ram_sv(ikl) =0. 134 rah_sv(ikl) =0. 135 Fh__sv(ikl) =0. 136 dFh_sv(ikl) =0. 137 Evp_sv(ikl) =0. 138 EvT_sv(ikl) =0. 139 LSdzsv(ikl) =0. 140 Tsrfsv(ikl) =0. 141 sEX_sv(ikl,:) =0. 142 zzsnsv(ikl,:) =0. 143 psi_sv(ikl,:) =0. 144 Khydsv(ikl,:) =0. 145 EExcsv(ikl) =0. 146 147 115 148 END DO 116 149 -
LMDZ6/trunk/libf/phylmd/inlandsis/inlandsis.F
r3792 r3900 1 subroutine INLANDSIS(SnoMod,BloMod,jjtime )1 subroutine INLANDSIS(SnoMod,BloMod,jjtime,debut) 2 2 3 3 USE dimphy … … 173 173 USE VARySV 174 174 USE VARtSV 175 USE surface_data, only: iflag_tsurf_inlandsis 176 175 USE surface_data, ONLY: is_ok_z0h_rn, 176 . is_ok_density_kotlyakov, 177 . prescribed_z0m_snow, 178 . iflag_z0m_snow, 179 . iflag_tsurf_inlandsis, 180 . iflag_temp_inlandsis, 181 . discret_xf, buf_sph_pol,buf_siz_pol 177 182 178 183 IMPLICIT NONE … … 180 185 logical SnoMod 181 186 logical BloMod 187 logical debut 182 188 integer jjtime 183 189 … … 213 219 integer IceMsk,IcIndx(klonv) ! Ice / No Ice Mask 214 220 integer SnoMsk ! Snow / No Snow Mask 215 216 221 real roSMin,roSMax,roSn_1,roSn_2,roSn_3 ! Fallen Snow Density (PAHAUT) 217 222 real Dendr1,Dendr2,Dendr3 ! Fallen Snow Dendric.(GIRAUD) 218 223 real Spher1,Spher2,Spher3,Spher4 ! Fallen Snow Spheric.(GIRAUD) 219 224 real Polair ! Polar Snow Switch 220 real PorSno, Por_BS,Salt_f,PorRef !225 real PorSno,Salt_f,PorRef ! 221 226 c #sw real PorVol,rWater ! 222 227 c #sw real rusNEW,rdzNEW,etaNEW ! … … 244 249 real Z0m_Sn,Z0m_90 ! Snow Surface Roughness Length 245 250 real SnoWat ! Snow Layer Switch 246 c #RNreal rstar,alors !247 c #RNreal rstar0,rstar1,rstar2 !251 real rstar,alors ! 252 real rstar0,rstar1,rstar2 ! 248 253 real SameOK ! 1. => Same Type of Grains 249 254 real G1same ! Averaged G1, same Grains … … 263 268 real Sph_av ! Averaged Grain Spher. 264 269 real Den_av ! Averaged Grain Dendr. 265 real DendOK ! 1. => Average is Dendr.266 270 real G1diff ! Averaged G1, diff. Grains 267 271 real G2diff ! Averaged G2, diff. Grains … … 277 281 real tt_c,vv_c ! Critical param. 278 282 real tt_tmp,vv_tmp,vv_virt ! Temporary variables 279 logical density_kotlyakov ! .true. if Kotlyakov 1961280 283 real e_prad,e1pRad,A_Rad0,absg_V,absgnI,exdRad ! variables for SoSosv calculations 281 284 real zm1, zm2, coefslope ! variables for surface temperature extrapolation 282 285 ! for Aeolian erosion and blowing snow 286 integer nit ,iit 287 real Fac ! Correc. factor for drift ratio 288 real dusuth,signus 289 real sss__F,sss__N 290 real sss__K,sss__G 291 real us_127,us_227,us_327,us_427,us_527 292 real VVa_OK, usuth0 293 real ssstar 294 real SblPom 295 real rCd10n ! Square root of drag coefficient 296 real DendOK ! Dendricity Switch 297 real SaltOK ! Saltation Switch 298 real MeltOK ! Saltation Switch (Melting Snow) 299 real SnowOK ! Pack Top Switch 300 real SaltM1,SaltM2,SaltMo,SaltMx ! Saltation Parameters 301 real ShearX, ShearS ! Arg. Max Shear Stress 302 real Por_BS ! Snow Porosity 303 real Salt_us ! New thresh.friction velocity u*t 304 real Fac_Mo,ArguSi,FacRho ! Numerical factors for u*t 305 real SaltSI(klonv,0:nsno) ! Snow Drift Index ! 306 real MIN_Mo ! Minimum Mobility Fresh Fallen * 307 character*3 qsalt_param ! Switch for saltation flux param. 308 character*3 usth_param ! Switch for u*t param 283 309 284 310 … … 287 313 288 314 data T__Min / 200.00/ ! Minimum realistic Temperature 289 data TaPole / 26 3.15/ ! Maximum Polar Temperature290 data roSMin / 30. / ! Minimum Snow Density315 data TaPole / 268.15/ ! Maximum Polar Temperature (value from C. Agosta) 316 data roSMin / 300. / ! Minimum Snow Density 291 317 data roSMax / 400. / ! Max Fresh Snow Density 292 318 data tt_c / -2.0 / ! Critical Temp. (degC) … … 305 331 data EmiWat / 0.99999999/ ! Emissivity of a Water Area 306 332 data EmiSno / 0.99999999/ ! Emissivity of Snow 333 307 334 308 335 ! DATA Emissivities ! Pielke, 1984, pp. 383,409 … … 321 348 data Z0_ICE/ 0.0010/ ! Sea-Ice Z0 = 0.0010 m (Andreas) 322 349 ! ! (Ice Station Weddel -- ISW) 350 ! for aerolian erosion 351 data SblPom/ 1.27/ ! Lower Boundary Height Parameter 352 C + ! for Suspension 353 C + ! Pommeroy, Gray and Landine 1993, 354 C + ! J. Hydrology, 144(8) p.169 355 data nit / 5 / ! us(is0,uth) recursivity: Nb Iterations 356 cc#AE data qsalt_param/"bin"/ ! saltation part. conc. from Bintanja 2001 (p 357 data qsalt_param/"pom"/ ! saltation part. conc. from Pomeroy and Gray 358 cc#AE data usth_param/"lis"/ ! u*t from Liston et al. 2007 359 data usth_param/"gal"/ ! u*t from Gallee et al. 2001 360 data SaltMx/-5.83e-2/ 361 323 362 vk2 = vonKrm * vonKrm ! Square of Von Karman Constant 324 363 … … 352 391 353 392 354 ! Blowing Particles Threshold Friction velocity 355 ! ============================================= 356 357 c #AE usthSV(ikl) = 1.0e+2 358 ! END DO 359 !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 360 361 362 363 364 ! Contribution of Snow to the Surface Snow Pack 365 ! ============================================= 366 367 IF (SnoMod) THEN 368 369 370 371 C +--Blowing Snow 372 C + ------------ 373 374 IF (BloMod) then 375 if (klonv.eq.1) then 393 394 395 396 IF (SnoMod) THEN 397 398 399 C +--Aeolian erosion and Blowing Snow 400 C +================================== 401 402 403 404 DO ikl=1,knonv 405 usthSV(ikl) = 1.0e+2 406 END DO 407 408 409 IF (BloMod) THEN 410 411 if (klonv.eq.1) then 376 412 if(isnoSV(1).ge.2 .and. 377 . TsisSV(1,max(1,isnoSV(1)))<273. .and.378 . ro__SV(1,max(1,isnoSV(1)))<500. .and.379 . eta_SV(1,max(1,isnoSV(1)))<epsi) then413 . TsisSV(1,max(1,isnoSV(1)))<273. .and. 414 . ro__SV(1,max(1,isnoSV(1)))<500. .and. 415 . eta_SV(1,max(1,isnoSV(1)))<epsi) then 380 416 C + ********** 381 417 call SISVAT_BSn … … 384 420 call SISVAT_BSn 385 421 C + ********** 386 endif 387 ENDIF 388 389 390 422 endif 423 424 425 426 427 428 429 430 ! Calculate threshold erosion velocity for next time step 431 ! Unlike in sisvat, computation is of threshold velocity made here (instead of sisvaesbl) 432 ! since we do not use sisvatesbl for the coupling with LMDZ 433 434 C +--Computation of threshold friction velocity for snow erosion 435 C --------------------------------------------------------------- 436 437 rCd10n = 1. / 26.5 ! Vt / u*t = 26.5 438 ! Budd et al. 1965, Antarct. Res. Series Fig.13 439 ! ratio developped during assumed neutral conditions 440 441 442 C +--Snow Properties 443 C + ~~~~~~~~~~~~~~~ 444 445 DO ikl = 1,knonv 446 447 isn = isnoSV(ikl) 448 449 450 451 DendOK = max(zero,sign(unun,epsi-G1snSV(ikl,isn) )) ! 452 SaltOK = min(1 , max(istdSV(2)-istoSV(ikl,isn),0)) ! 453 MeltOK = (unun ! 454 . -max(zero,sign(unun,TfSnow-epsi ! 455 . -TsisSV(ikl,isn) ))) ! Melting Snow 456 . * min(unun,DendOK ! 457 . +(1.-DendOK) ! 458 . *sign(unun, G2snSV(ikl,isn)-1.0)) ! 1.0 for 1mm 459 SnowOK = min(1 , max(isnoSV(ikl) +1 -isn ,0)) ! Snow Switch 460 461 G1snSV(ikl,isn) = SnowOK * G1snSV(ikl,isn) 462 . + (1.- SnowOK)*min(G1snSV(ikl,isn),G1_dSV) 463 G2snSV(ikl,isn) = SnowOK * G2snSV(ikl,isn) 464 . + (1.- SnowOK)*min(G2snSV(ikl,isn),G1_dSV) 465 466 SaltOK = min(unun, SaltOK + MeltOK) * SnowOK 467 468 469 C +--Mobility Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.) 470 C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 471 SaltM1 = -0.750e-2 * G1snSV(ikl,isn) 472 . -0.500e-2 * G2snSV(ikl,isn)+ 0.500e00 !dendritic case 473 C + CAUTION: Guyomarc'h & Merindol Dendricity Sign is + 474 C + ^^^^^^^^ MAR Dendricity Sign is - 475 SaltM2 = -0.833d-2 * G1snSV(ikl,isn) 476 . -0.583d-2 * G2snSV(ikl,isn)+ 0.833d00 !non-dendritic case 477 478 c SaltMo = (DendOK * SaltM1 + (1.-DendOK) * SaltM2 ) 479 SaltMo = 0.625 !SaltMo pour d=s=0.5 480 481 !weighting SaltMo with surface snow density (Vionnet et al. 2012) 482 cc#AE FacRho = 1.25 - 0.0042 * ro__SV(ikl,isn) 483 cc#AE SaltMo = 0.34 * SaltMo + 0.66 * FacRho !needed for polar snow 484 MIN_Mo = 0. 485 c SaltMo = max(SaltMo,MIN_Mo) 486 c SaltMo = SaltOK * SaltMo + (1.-SaltOK) * min(SaltMo,SaltMx) 487 c #TUNE SaltMo = SaltOK * SaltMo - (1.-SaltOK) * 0.9500 488 SaltMo = max(SaltMo,epsi-unun) 489 490 C +--Influence of Density on Threshold Shear Stress 491 C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 492 Por_BS = 1. - 300. / ro_Ice 493 ShearS = Por_BS / (1.-Por_BS) 494 C +... SheaBS = Arg(sqrt(shear = max shear stress in snow)): 495 C + shear = 3.420d00 * exp(-(Por_BS +Por_BS) 496 C + . /(unun -Por_BS)) 497 C + SheaBS : see de Montmollin (1978), 498 C + These Univ. Sci. Medic. Grenoble, Fig. 1 p. 124 499 500 C +--Snow Drift Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.) 501 C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 502 ArguSi = -0.085 *us__SV(ikl)/rCd10n 503 !V=u*/sqrt(CD) eqs 2 to 4 Gallee et al. 2001 504 505 SaltSI(ikl,isn) = -2.868 * exp(ArguSi) + 1 + SaltMo 506 507 508 C +--Threshold Friction Velocity 509 C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 510 if(ro__SV(ikl,isn)>300.) then 511 Por_BS = 1.000 - ro__SV(ikl,isn) /ro_Ice 512 else 513 Por_BS = 1.000 - 300. /ro_Ice 514 endif 515 516 ShearX = Por_BS/max(epsi,1.-Por_BS) 517 Fac_Mo = exp(-ShearX+ShearS) 518 C + Gallee et al., 2001 eq 5, p5 519 520 if (usth_param .eq. "gal") then 521 Salt_us = (log(2.868) - log(1 + SaltMo)) * rCd10n/0.085 522 Salt_us = Salt_us * Fac_Mo 523 C +... Salt_us : Extension of Guyomarc'h & Merindol 1998 with 524 C +... de Montmollin (1978). Gallee et al. 2001 525 endif 526 527 if (usth_param .eq. "lis") then !Liston et al. 2007 528 if(ro__SV(ikl,isn)>300.) then 529 Salt_us = 0.005*exp(0.013*ro__SV(ikl,isn)) 530 else 531 Salt_us = 0.01*exp(0.003*ro__SV(ikl,isn)) 532 endif 533 endif 534 535 SnowOK = 1 -min(1,iabs(isn-isnoSV(ikl))) !Switch new vs old snow 536 537 usthSV(ikl) = SnowOK * (Salt_us) 538 . + (1.-SnowOK)* usthSV(ikl) 539 540 END DO 541 542 543 544 ! Feeback between blowing snow turbulent Scale u* (commented here 545 ! since ustar is an input variable (not in/out) of inlandsis) 546 ! ----------------------------------------------------------------- 547 548 549 ! VVa_OK = max(0.000001, VVaSBL(ikl)) 550 ! sss__N = vonkar * VVa_OK 551 ! sss__F = (sqrCm0(ikl) - psim_z + psim_0) 552 ! usuth0 = sss__N /sss__F ! u* if NO Blow. Snow 553 554 ! sss__G = 0.27417 * gravit 555 556 ! ! ______________ _____ 557 ! ! Newton-Raphson (! Iteration, BEGIN) 558 ! ! ~~~~~~~~~~~~~~ ~~~~~ 559 ! DO iit=1,nit 560 ! sss__K = gravit * r_Turb * A_Turb *za__SV(ikl) 561 ! . *rCDmSV(ikl)*rCDmSV(ikl) 562 ! . /(1.+0.608*QaT_SV(ikl)-qsnoSV(ikl)) 563 ! us_127 = exp( SblPom *log(us__SV(ikl))) 564 ! us_227 = us_127 * us__SV(ikl) 565 ! us_327 = us_227 * us__SV(ikl) 566 ! us_427 = us_327 * us__SV(ikl) 567 ! us_527 = us_427 * us__SV(ikl) 568 569 ! us__SV(ikl) = us__SV(ikl) 570 ! . - ( us_527 *sss__F /sss__N 571 ! . - us_427 572 ! . - us_227 *qsnoSV(ikl)*sss__K 573 ! . + (us__SV(ikl)*us__SV(ikl)-usthSV(ikl)*usthSV(ikl))/sss__G) 574 ! . /( us_427*5.27*sss__F /sss__N 575 ! . - us_327*4.27 576 ! . - us_127*2.27*qsnoSV(ikl)*sss__K 577 ! . + us__SV(ikl)*2.0 /sss__G) 578 579 ! us__SV(ikl)= min(us__SV(ikl),usuth0) 580 ! us__SV(ikl)= max(us__SV(ikl),epsi ) 581 ! rCDmSV(ikl)= us__SV(ikl)/VVa_OK 582 ! ! #AE sss__F = vonkar /rCDmSV(ikl) 583 ! ENDDO 584 585 ! ! ______________ ___ 586 ! ! Newton-Raphson (! Iteration, END ) 587 ! ! ~~~~~~~~~~~~~~ ~~~ 588 589 ! us_127 = exp( SblPom *log(us__SV(ikl))) 590 ! us_227 = us_127 * us__SV(ikl) 591 592 ! ! Momentum Turbulent Scale u*: 0-Limit in case of no Blow. Snow 593 ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 594 ! dusuth = us__SV(ikl) - usthSV(ikl) ! u* - uth* 595 ! signus = max(sign(unun,dusuth),zero) ! 1 <=> u* - uth* > 0 596 ! us__SV(ikl) = ! 597 ! . us__SV(ikl) *signus + ! u* (_BS) 598 ! . usuth0 ! u* (nBS) 599 ! . *(1.-signus) ! 600 601 602 603 604 ! Blowing Snow Turbulent Scale ss* 605 ! --------------------------------------- 606 607 hSalSV(ikl) = 8.436e-2 * us__SV(ikl)**SblPom 608 609 if (qsalt_param .eq. "pom") then 610 qSalSV(ikl) = (us__SV(ikl)**2 - usthSV(ikl)**2) *signus 611 . / (hSalSV(ikl) * gravit * us__SV(ikl) * 3.25) 612 endif 613 614 if (qsalt_param .eq. "bin") then 615 qSalSV(ikl) = (us__SV(ikl) * us__SV(ikl) 616 . -usthSV(ikl) * usthSV(ikl))*signus 617 . * 0.535 / (hSalSV(ikl) * gravit) 618 endif 619 620 qSalSV(ikl) = qSalSV(ikl)/rht_SV(ikl) ! conversion kg/m3 to kg/kg 621 622 ssstar = rCDmSV(ikl) * (qsnoSV(ikl) - qSalSV(ikl)) 623 . * r_Turb !Bintanja 2000, BLM 624 !r_Turb compensates for an overestim. of the blown snow part. fall velocity 625 626 uss_SV(ikl) = min(zero , us__SV(ikl) *ssstar) 627 uss_SV(ikl) = max(-0.0001 , uss_SV(ikl)) 628 629 630 631 632 ENDIF ! BloMod 633 634 C + ------------------------------------------------------ 391 635 C +--Buffer Layer 392 C + ------------ 636 C + ----------------------------------------------------- 393 637 394 638 DO ikl=1,knonv … … 414 658 c #NP. 104. *sqrt( max( VV10SV(ikl)-6.0,0.0))) ! Kotlyakov (1961) 415 659 416 density_kotlyakov = .true.417 c #AC density_kotlyakov = .false. !C.Agosta snow densisty as if BS is on b 660 ! C.Agosta option for snow density, same as for BS i.e. 661 ! is_ok_density_kotlyakov=.false. 418 662 c #BS density_kotlyakov = .false. !C.Amory BS 2018 419 663 C + ... Fallen Snow Density, Adapted for Antarctica 420 if ( density_kotlyakov) then664 if (is_ok_density_kotlyakov) then 421 665 tt_tmp = TaT_SV(ikl)-TfSnow 422 666 !vv_tmp = VV10SV(ikl) … … 452 696 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 453 697 454 c #BS Bros_N = frsno 455 c #BS ro_new = ro__SV(ikl,max(1,isnoSV(ikl))) 456 c #BS ro_new = max(Bros_N,min(roBdSV,ro_new)) 457 c #BS Fac = 1-((ro__SV(ikl,max(1,isnoSV(ikl))) 458 c #BS. -roBdSV)/(500.-roBdSV)) 459 c #BS Fac = max(0.,min(1.,Fac)) 460 c #BS dsnbSV(ikl) = Fac*dsnbSV(ikl) 461 c #BS Bros_N = Bros_N * (1.0-dsnbSV(ikl)) 462 c #BS. + ro_new * dsnbSV(ikl) 463 698 if (BloMod) then 699 Bros_N = frsno 700 ro_new = ro__SV(ikl,max(1,isnoSV(ikl))) 701 ro_new = max(Bros_N,min(roBdSV,ro_new)) 702 Fac = 1-((ro__SV(ikl,max(1,isnoSV(ikl))) 703 . -roBdSV)/(500.-roBdSV)) 704 Fac = max(0.,min(1.,Fac)) 705 dsnbSV(ikl) = Fac*dsnbSV(ikl) 706 Bros_N = Bros_N * (1.0-dsnbSV(ikl)) 707 . + ro_new * dsnbSV(ikl) 708 endif 464 709 465 710 … … 480 725 . max(Spher1*VV__SV(ikl)+Spher2, ! Sphericity 481 726 . Spher3 )) ! 727 ! EV: now control buf_sph_pol and bug_siz_pol in physiq.def 482 728 Buf_G1 = (1. - Polair) * Buf_G1 ! Temperate Snow 483 . + Polair * G1_dSV ! Polar Snow729 . + Polair * buf_sph_pol ! Polar Snow 484 730 Buf_G2 = (1. - Polair) * Buf_G2 ! Temperate Snow 485 . + Polair * ADSdSV ! PolarSnow731 . + Polair * buf_siz_pol ! Polar Snow 486 732 G1 = Buf_G1 ! NO Blown Snow 487 733 G2 = Buf_G2 ! NO Blown Snow 488 734 489 735 736 737 IF (BloMod) THEN 738 490 739 ! S.1. Meme Type de Neige / same Grain Type 491 740 ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 492 c #BS SameOK = max(zero, 493 c #BS. sign(unun, Buf_G1 *G1_dSV 494 c #BS. - eps_21 )) 495 c #BS G1same = ((1.0-dsnbSV(ikl))*Buf_G1+dsnbSV(ikl) *G1_dSV) 496 c #BS G2same = ((1.0-dsnbSV(ikl))*Buf_G2+dsnbSV(ikl) *ADSdSV) 741 742 SameOK = max(zero, 743 . sign(unun, Buf_G1 *G1_dSV 744 . - eps_21 )) 745 G1same = ((1.0-dsnbSV(ikl))*Buf_G1+dsnbSV(ikl) *G1_dSV) 746 G2same = ((1.0-dsnbSV(ikl))*Buf_G2+dsnbSV(ikl) *ADSdSV) 497 747 ! Blowing Snow Properties: G1_dSV, ADSdSV 498 748 499 749 ! S.2. Types differents / differents Types 500 750 ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 501 c #BStyp__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic502 c #BSzroNEW = typ__1 *(1.0-dsnbSV(ikl)) ! fract.Dendr.Lay.503 c #BS.+ (1.-typ__1) * dsnbSV(ikl) !504 c #BSG1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay.505 c #BS.+ (1.-typ__1) *G1_dSV !506 c #BSG2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay.507 c #BS.+ (1.-typ__1) *ADSdSV !508 c #BSzroOLD = (1.-typ__1) *(1.0-dsnbSV(ikl)) ! fract.Spher.Lay.509 c #BS.+ typ__1 * dsnbSV(ikl) !510 c #BSG1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay.511 c #BS.+ typ__1 *G1_dSV !512 c #BSG2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay.513 c #BS.+ typ__1 *ADSdSV !514 c #BSSizNEW = -G1_NEW *DDcdSV/G1_dSV ! Size Dendr.Lay.515 c #BS. +(1.+G1_NEW /G1_dSV)!516 c #BS. *(G2_NEW *DScdSV/G1_dSV!517 c #BS. +(1.-G2_NEW /G1_dSV)*DFcdSV)!518 c #BSSphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay.519 c #BSSizOLD = G2_OLD ! Size Spher.Lay.520 c #BSSphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay.521 c #BS Siz_av = (zroNEW*SizNEW+zroOLD*SizOLD)! Averaged Size522 c #BS Sph_av = min( zroNEW*SphNEW+zroOLD*SphOLD!523 c #BS., unun) ! Averaged Sphericity524 c #BS Den_av = min((Siz_av -( Sph_av *DScdSV!525 c #BS. +(1.-Sph_av)*DFcdSV))!526 c #BS. / (DDcdSV -( Sph_av *DScdSV!527 c #BS. +(1.-Sph_av)*DFcdSV))!528 c #BS. , unun)!529 c #BSDendOK = max(zero, !530 c #BS. sign(unun, Sph_av *DScdSV ! Small Grains531 c #BS. +(1.-Sph_av)*DFcdSV ! Faceted Grains532 c #BS. - Siz_av )) !751 typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic 752 zroNEW = typ__1 *(1.0-dsnbSV(ikl)) ! fract.Dendr.Lay. 753 . + (1.-typ__1) * dsnbSV(ikl) ! 754 G1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay. 755 . + (1.-typ__1) *G1_dSV ! 756 G2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay. 757 . + (1.-typ__1) *ADSdSV ! 758 zroOLD = (1.-typ__1) *(1.0-dsnbSV(ikl)) ! fract.Spher.Lay. 759 . + typ__1 * dsnbSV(ikl) ! 760 G1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay. 761 . + typ__1 *G1_dSV ! 762 G2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay. 763 . + typ__1 *ADSdSV ! 764 SizNEW = -G1_NEW *DDcdSV/G1_dSV ! Size Dendr.Lay. 765 . +(1.+G1_NEW /G1_dSV) ! 766 . *(G2_NEW *DScdSV/G1_dSV ! 767 . +(1.-G2_NEW /G1_dSV)*DFcdSV) ! 768 SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay. 769 SizOLD = G2_OLD ! Size Spher.Lay. 770 SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay. 771 Siz_av = (zroNEW*SizNEW+zroOLD*SizOLD) ! Averaged Size 772 Sph_av = min( zroNEW*SphNEW+zroOLD*SphOLD ! 773 . , unun) ! Averaged Sphericity 774 Den_av = min((Siz_av -( Sph_av *DScdSV ! 775 . +(1.-Sph_av)*DFcdSV)) ! 776 . / (DDcdSV -( Sph_av *DScdSV ! 777 . +(1.-Sph_av)*DFcdSV)) ! 778 . , unun) ! 779 DendOK = max(zero, ! 780 . sign(unun, Sph_av *DScdSV ! Small Grains 781 . +(1.-Sph_av)*DFcdSV ! Faceted Grains 782 . - Siz_av )) ! 533 783 C +... REMARQUE: le type moyen (dendritique ou non) depend 534 784 C + ^^^^^^^^ de la comparaison avec le diametre optique … … 538 788 C + of a recent snow having zero dendricity 539 789 540 c #BS G1diff =( -DendOK *Den_av 541 c #BS. +(1.-DendOK)*Sph_av) *G1_dSV 542 c #BS G2diff = DendOK *Sph_av *G1_dSV 543 c #BS. +(1.-DendOK)*Siz_av 544 c #BS G1 = SameOK *G1same 545 c #BS. +(1.-SameOK)*G1diff 546 c #BS G2 = SameOK *G2same 547 c #BS. +(1.-SameOK)*G2diff 548 790 G1diff =( -DendOK *Den_av 791 . +(1.-DendOK)*Sph_av) *G1_dSV 792 G2diff = DendOK *Sph_av *G1_dSV 793 . +(1.-DendOK)*Siz_av 794 G1 = SameOK *G1same 795 . +(1.-SameOK)*G1diff 796 G2 = SameOK *G2same 797 . +(1.-SameOK)*G2diff 798 ENDIF 799 549 800 550 801 … … 634 885 . /max(epsi,BrosSV(ikl))!& [m w.e.] -> [m] 635 886 636 637 887 638 888 END DO … … 640 890 641 891 642 ! Snow Pack Discretization 643 ! ======================== 644 645 ! ********** 892 ! Snow Pack Discretization(option XF in MAR) 893 ! ========================================== 894 895 896 if (discret_xf.AND.klonv.eq.1) then 897 898 if(isnoSV(1).ge.1.or.NLaysv(1).ge.1) then 899 C + ********** 900 call SISVAT_zSn 901 C + ********** 902 endif 903 else 904 C + ********** 646 905 call SISVAT_zSn 647 ! ********** 648 649 ! ********** 906 C + ********** 907 endif 908 909 C + ********** 650 910 ! #ve call SISVAT_wEq('_zSn ',0) 651 ! ********** 652 653 911 C + ********** 654 912 655 913 ! Add a new Snow Layer … … 664 922 TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl)) 665 923 . + min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl) 666 667 924 ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl)) 668 925 . + Brossv(ikl) * NLaysv(ikl) … … 699 956 700 957 701 END IF 958 END IF ! SnoMod 702 959 703 960 … … 740 997 ! ============================= 741 998 !Etienne: as in inlandis we do not call vgopt, we need to define 742 !the albedo 999 !the albedo alb_SV and to calculate the 743 1000 !absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv 744 1001 … … 810 1067 811 1068 812 ! Aerodynamic Resistance 813 ! ^^^^^^^^^^^^^^^^^^^^^^ 814 815 816 DO ikl=1,knonv 817 ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6)) 818 rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6)) 819 END DO 1069 ! Aerodynamic Resistance (calculated from drags given by LMDZ) 1070 ! Commented because already calculated in surf_inlandsis_mod 1071 ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1072 ! DO ikl=1,knonv 1073 ! ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6)) 1074 ! rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6)) 1075 ! END DO 820 1076 821 1077 … … 825 1081 826 1082 827 if (iflag_t surf_inlandsis .eq. 0) then1083 if (iflag_temp_inlandsis .eq. 0) then 828 1084 829 1085 call SISVAT_TSo 830 1086 831 1087 else 1088 DO ikl=1,knonv 1089 Tsf_SV(ikl)=Tsrfsv(ikl) 1090 END DO 832 1091 833 1092 call SISVAT_TS2 … … 938 1197 ! Surface Temperature 939 1198 ! ^^^^^^^^^^^^^^^^^^^^ 940 ! Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl)) 941 1199 1200 IF (iflag_tsurf_inlandsis .EQ. 0) THEN 1201 1202 Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl)) 1203 1204 ELSE IF (iflag_tsurf_inlandsis .GT. 0) THEN 942 1205 ! Etienne: extrapolation from the two uppermost levels: 943 1206 … … 959 1222 960 1223 961 END DO 962 1224 ELSE !(default) 1225 1226 Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl)) 1227 1228 END IF 1229 1230 1231 END DO 963 1232 964 1233 ! Snow Pack Properties (sphericity, dendricity, size) … … 967 1236 IF (SnoMod) THEN 968 1237 969 ! ********** 1238 if (discret_xf .AND. klonv.eq.1) then 1239 if(isnoSV(1).ge.1) then 1240 C + ********** 1241 call SISVAT_GSn 1242 C + ********** 1243 endif 1244 else 1245 C + ********** 970 1246 call SISVAT_GSn 971 ! ********** 972 973 ! ********** 974 ! #ve call SISVAT_wEq('_GSn ',0) 975 ! ********** 976 1247 C + ********** 1248 endif 977 1249 978 1250 … … 990 1262 C +--Roughness Length for Momentum 991 1263 C + ----------------------------- 1264 1265 ! ETIENNE WARNING: changes have been made wrt original SISVAT 992 1266 993 1267 C +--Land+Sea-Ice / Ice-free Sea Mask 994 1268 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 995 DO ikl=1,k lonv1269 DO ikl=1,knonv 996 1270 IcIndx(ikl) = 0 997 1271 ENDDO 998 1272 DO isn=1,nsno 999 DO ikl=1,klonv 1273 DO ikl=1,knonv 1274 1000 1275 IcIndx(ikl) = max(IcIndx(ikl), 1001 . 1002 . 1003 . 1276 . isn*max(0, 1277 . sign(1, 1278 . int(ro__SV(ikl,isn)-900.)))) 1004 1279 ENDDO 1005 1280 ENDDO 1006 1281 1007 DO ikl=1,k lonv1282 DO ikl=1,knonv 1008 1283 LISmsk = 1. ! in inlandsis, land only 1009 1284 IceMsk = max(0,sign(1 ,IcIndx(ikl)-1) ) 1010 1285 SnoMsk = max(min(isnoSV(ikl)-iiceSV(ikl),1),0) 1011 1286 1012 1013 1014 Z0mLnd =max( Z0_ICE , 5.e-5 ) ! Min set := Z0 on *1015 1287 1016 1288 C +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8) 1017 1289 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ 1018 1290 Z0m_nu = 5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03)) 1019 1291 1020 1292 C +--Z0 Saltat.Regime over Snow (Gallee et al., 2001, BLM 99 (19) p.11) 1021 1293 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ 1294 1022 1295 u2star = us__SV(ikl) *us__SV(ikl) 1023 1296 Z0mBSn = u2star *0.536e-3 - 61.8e-6 1024 1297 Z0mBSn = max(Z0mBS0 ,Z0mBSn) 1025 1298 1026 1299 C +--Z0 Smooth + Saltat. Regime 1027 1300 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ 1028 1301 Z0enSV(ikl) = Z0m_nu 1029 1302 . + Z0mBSn 1030 1031 C +--Rough Snow Surface Roughness Length (Typical Value) 1032 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1033 c #tz Z0m_Sn = 0.250e-3 ! Andreas 1995, CRREL Report 95-16, fig.1&p.2 1034 ! z0r~(10-d)*exp(-vonkar/sqrt(1.5e-03))-5.e-5 1035 Z0m_Sn = 2.000e-3 ! Calibration of MAR 1036 c #TZ Z0m_Sn = 1.000e-3 ! Exemple Tuning in RACMO 1037 c #TZ Z0m_Sn = 0.500e-3 ! Exemple Tuning in MAR 1038 1303 1304 1305 ! Calculation of snow roughness length 1306 !===================================== 1307 IF (iflag_z0m_snow .EQ. 0) THEN 1308 1309 Z0m_Sn=prescribed_z0m_snow 1310 1311 ELSE IF (iflag_z0m_snow .EQ. 1) THEN 1312 1313 Z0m_Sn=Z0enSV(ikl) 1314 1315 ELSE IF (iflag_z0m_snow .EQ. 2) THEN 1316 1039 1317 C +--Rough Snow Surface Roughness Length (Variable Sastrugi Height) 1040 1318 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ … … 1045 1323 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1046 1324 ! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013 1325 1326 1047 1327 coefa = 0.1658 !0.1862 !Ant 1048 1328 coefb = -50.3869 !-55.7718 !Ant … … 1055 1335 coefc = log(z03/z02)/(ta3-ta2) 1056 1336 coefd = log(z03)-coefc*ta3 1337 1057 1338 if (TaT_SV(ikl) .lt. ta1) then 1058 1339 Z0_obs = z01 … … 1066 1347 endif 1067 1348 1068 1069 ! pour le moment, on choisit une valeur fixe 1070 Z0_obs = 1.000e-3 1071 1072 cCA Snow roughness lenght deduced from observations 1073 cCA (parametrization if no Blowing Snow) 1074 cCA ----------------------------------- C. Agosta 09-2016 ----- 1075 cCA Substract Z0enSV(ikl) because re-added later in Z0mnSV(ikl) 1076 Z0m_Sn = Z0_obs - Z0enSV(ikl) 1077 cCA ----------------------------------------------------------- 1078 1079 param = Z0_obs/1. ! param(s) | 1.(m/s)=TUNING 1080 1349 Z0m_Sn=Z0_obs 1350 1351 1352 ELSE 1353 1354 Z0m_Sn=0.500e-3 ! default=0.500e-3m (tuning of MAR) 1355 1356 ENDIF 1357 1358 1359 1360 ! param = Z0_obs/1. ! param(s) | 1.(m/s)=TUNING 1081 1361 c #SZ Z0Sa_N = (us__SV(ikl) -0.2)*param ! 0.0001=TUNING 1082 1362 c #SZ. * max(zero,sign(unun,TfSnow-eps9 … … 1109 1389 c #ZN Z0enSV(ikl) = max(Z0enSV(ikl), Z0m_nu) 1110 1390 1391 1111 1392 C +--Z0 Smooth Regime over Snow (Andreas etAl., 2004 1112 1393 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf) … … 1132 1413 c #ZA Z0m_Sn = DDs_SV(ikl)* Z0m_90 / 45. 1133 1414 c #ZA. - DDs_SV(ikl)*DDs_SV(ikl)* Z0m_90 /(90.*90.) 1134 1135 C +--Z0 (Erosion) over Snow (instantaneous or time average) 1415 1416 1417 1418 1419 C +--Z0 (Erosion) over Snow (instantaneous) 1136 1420 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ 1137 1421 Z0e_SV(ikl) = Z0enSV(ikl) 1138 Z0e_SV(ikl) = Z0emSV(ikl) 1139 1140 C +--Momentum Roughness Length 1141 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Contribution of 1142 Z0mnSV(ikl) = Z0mLnd ! land Form 1143 . + (Z0m_Sn ! Sastrugi Form 1144 . + Z0enSV(ikl)) *SnoMsk ! Snow Erosion 1422 1423 C +--Momentum Roughness Length (Etienne: changes wrt original SISVAT) 1424 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1425 Z0mnSV(ikl) = Z0m_nu *(1-SnoMsk) ! Ice z0 1426 . + (Z0m_Sn)*SnoMsk ! Snow Sastrugi Form and Snow Erosion 1145 1427 1146 1428 … … 1154 1436 c #GL. /(920.00 -600.))) ! 1155 1437 1156 C +--Mom. Roughness Length, Instantaneous OR Box Moving Average in Time1157 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^1438 C +--Mom. Roughness Length, Instantaneous 1439 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1158 1440 Z0m_SV(ikl) = Z0mnSV(ikl) ! Z0mnSV instant. 1159 ! Z0m_SV(ikl) = Z0mmSV(ikl) ! Z0mnSV Average1160 1161 C +--Corrected Threshold Friction Velocity before Erosion ! Marticorena and1162 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Bergametti 19951163 ! not used anymore since Marticorena and Bergametti disabled !CK 18/07/20181164 cc #BS Z0e_SV(ikl) = min(Z0m_SV(ikl),Z0e_SV(ikl)) !1165 cc #MB f_eff= log(0.35*(0.1 /Z0e_SV(ikl))**0.8) ! JGR 1001166 cc #MB f_eff=1.-(log( Z0m_SV(ikl)/Z0e_SV(ikl) ))! (20) p. 164201167 cc #MB. /(max( f_eff ,epsi ))! p.16426 2nd ?1168 cc #MB f_eff= max( f_eff ,epsi )! CONTROL1169 cc #MB f_eff=1.0 -(1.0 - f_eff) /5.00 ! TUNING1170 cc #MB f_eff= min( f_eff ,1.00 )!1171 cc #MB usthSV(ikl) = usthSV(ikl)/f_eff !1172 1173 1441 1174 1442 … … 1177 1445 1178 1446 Z0hnSV(ikl) = Z0mnSV(ikl)/ 7.4 1179 c #SH Z0hnSV(ikl) = Z0mnSV(ikl)/100.0 1180 C + Z0h = Z0m /100.0 over the Sahel 1181 C + (Taylor & Clark, QJRMS 127,p864) 1182 1183 c #RN rstar = Z0mnSV(ikl) * us__SV(ikl) / akmol 1184 c #RN rstar = max(epsi,min(rstar,thous)) 1185 c #RN alors = log(rstar) 1186 c #RN rstar0 = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar)) 1187 c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar))) 1188 c #RN. *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar)) 1189 c #RN. + 0.317e0 1190 c #RN. *(1. - max(zero,sign(unun,2.500e0 - rstar)))) 1191 c #RN rstar1 = 0. * max(zero,sign(unun,0.135e0 - rstar)) 1192 c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar))) 1193 c #RN. *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar)) 1194 c #RN. - 0.565 1195 c #RN. *(1. - max(zero,sign(unun,2.500e0 - rstar)))) 1196 c #RN rstar2 = 0. * max(zero,sign(unun,0.135e0 - rstar)) 1197 c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar))) 1198 c #RN. *(0. * max(zero,sign(unun,2.500e0 - rstar)) 1199 c #RN. - 0.183 1200 c #RN. *(unun - max(zero,sign(unun,2.500e0 - rstar)))) 1201 1202 cXF #RN does not work over bare ice 1203 cXF MAR is then too warm and not enough melt 1204 1205 c #RN if(ro__SV(ikl,isnoSV(ikl))>50 1206 c #RN. .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then 1207 1208 c #RN Z0hnSV(ikl) = max(zero 1209 c #RN. , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)) 1210 c #RN. * exp(rstar0+rstar1*alors+rstar2*alors*alors) 1211 c #RN. * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero 1212 c #RN. , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))) 1213 1214 c #RN endif 1447 1448 IF (is_ok_z0h_rn) THEN 1449 1450 rstar = Z0mnSV(ikl) * us__SV(ikl) / akmol 1451 rstar = max(epsi,min(rstar,R_1000)) 1452 alors = log(rstar) 1453 rstar0 = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar)) 1454 . +(1. - max(zero,sign(unun,0.135e0 - rstar))) 1455 . *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar)) 1456 . + 0.317e0 1457 . *(1. - max(zero,sign(unun,2.500e0 - rstar)))) 1458 rstar1 = 0. * max(zero,sign(unun,0.135e0 - rstar)) 1459 . +(1. - max(zero,sign(unun,0.135e0 - rstar))) 1460 . *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar)) 1461 . - 0.565 1462 . *(1. - max(zero,sign(unun,2.500e0 - rstar)))) 1463 rstar2 = 0. * max(zero,sign(unun,0.135e0 - rstar)) 1464 . +(1. - max(zero,sign(unun,0.135e0 - rstar))) 1465 . *(0. * max(zero,sign(unun,2.500e0 - rstar)) 1466 . - 0.183 1467 . *(unun - max(zero,sign(unun,2.500e0 - rstar)))) 1468 1469 1470 1471 !XF #RN (is_ok_z0h_rn) does not work well over bare ice 1472 !XF MAR is then too warm and not enough melt 1473 1474 if(ro__SV(ikl,isnoSV(ikl))>50 1475 . .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then 1476 1477 Z0hnSV(ikl) = max(zero 1478 . , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)) 1479 . * exp(rstar0+rstar1*alors+rstar2*alors*alors) 1480 . * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero 1481 . , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))) 1482 1483 endif 1484 1485 1486 ENDIF 1215 1487 1216 1488 Z0h_SV(ikl) = Z0hnSV(ikl) 1217 ! Z0h_SV(ikl) = Z0hmSV(ikl)1218 1489 1219 1490 -
LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_bsn.F
r3792 r3900 9 9 C | | 10 10 C | SISVAT_bsn computes the snow erosion mass according to both the | 11 C | theoretical maximum erosion amount computed in SISVATesbl and the|11 C | theoretical maximum erosion amount computed in inlandsis and the | 12 12 C | availability of snow (currently in the uppermost snow layer only) | 13 13 C | | 14 C | Preprocessing Option: SISVAT IO (not always a standard preprocess.) |15 C | ^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^ |16 C | FILE | CONTENT |17 C | ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |18 C | # stdout | #sb: OUTPUT of Snow Erosion |19 C | | unit 6, SubRoutine SISVAT_BSn **ONLY** |20 14 C +------------------------------------------------------------------------+ 21 15 -
LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_qsn.F
r3792 r3900 61 61 use VARxSV 62 62 use VARySV 63 use surface_data, only: is_ok_slush,opt_runoff_ac 64 63 65 64 66 IMPLICIT NONE … … 235 237 236 238 DO ikl=1,knonv 237 DO isn=min(nsno,isnoSV(ikl)+1),1,-1 239 240 DO isn=min(nsno,isnoSV(ikl)+1),1,-1 238 241 ! EV DO isn=nsno,1,-1 239 242 C +--Energy, store Previous Content … … 243 246 . + ro__SV(ikl,isn) * Cn_dSV * dTSnow 244 247 . * dzsnSV(ikl,isn) 245 246 Tsave = TsisSV(ikl,isn)247 248 248 TsisSV(ikl,isn) = TfSnow 249 249 … … 312 312 rdzNEW = WaFrez + rdzsno 313 313 ro__SV(ikl,isn) = rdzNEW /max(epsi, dzsnSV(ikl,isn)) 314 315 ! EV: condition on Enfrez316 ! if (EnFrez .eq. 0.) then317 318 TsisSV(ikl,isn) = Tsave319 ! else320 314 TsisSV(ikl,isn) = TfSnow 321 315 . + EnFrez /(Cn_dSV *max(epsi, rdzNEW) ) 322 ! end if323 316 EExcsv(ikl) = EExcsv(ikl) - EnFrez 324 317 wer_SV(ikl) = WaFrez … … 499 492 rusnew = rusnSV(ikl) * SWf_SV(ikl) 500 493 501 if(isnoSV(ikl)<=1 ) rusnew = 0.494 if(isnoSV(ikl)<=1 .OR. opt_runoff_ac) rusnew = 0. 502 495 !if(ivgtSV(ikl)>=1) rusnew = 0. 503 496 504 497 c #EU rusnew = 0. 505 c #AC rusnew = 0. 498 c #AC rusnew = 0. 499 506 500 RnofSV(ikl) = RnofSV(ikl) 507 501 . +(rusnSV(ikl) - rusnew ) / dt__SV … … 545 539 ENDDO 546 540 547 C +--Slush Formation ( CAUTION: ADD RunOff Possibility before Activation)541 C +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation) 548 542 C + --------------- ^^^^^^^ ^^^ 549 543 550 551 c #SU DO ikl=1,knonv 552 c #SU DO isn=1,isnoSV(ikl) 553 c #SU kSlush = min(1,max(0,isn+1-ispiSV(ikl))) ! Slush Switch 544 IF (is_ok_slush) THEN 545 546 DO ikl=1,knonv 547 DO isn=1,isnoSV(ikl) 548 kSlush = min(1,max(0,isn+1-ispiSV(ikl))) ! Slush Switch 554 549 555 550 C +--Available Additional Pore Volume [-] 556 551 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 557 c #SUPorVol = 1. - ro__SV(ikl,isn) ! [--]558 c #SU. *(1. - eta_SV(ikl,isn))/ ro_Ice !559 c #SU. - eta_SV(ikl,isn) !560 c #SU. *ro__SV(ikl,isn) / ro_Wat !561 c #SUPorVol = max(PorVol , zero ) !562 c #SUzWater = dzsnSV(ikl,isn) * PorVol * 1000. ! [mm] OR [kg/m2]563 c #SU. * (1. -SWS_SV(ikl) ! 0 <=> freezing564 c #SU. *(1 -min(1,iabs(isn-isnoSV(ikl))))) ! 1 <=> isn=isnoSV565 c #SUzSlush = min(rusnSV(ikl) , zWater) ! [mm] OR [kg/m2]566 c #SUro_new =(dzsnSV(ikl,isn) * ro__SV(ikl,isn) !567 c #SU. +zSlush ) !568 c #SU. / max(dzsnSV(ikl,isn) , epsi ) !569 c #SUif(ro_new<ro_Ice+20) then ! MAX 940kg/m3 !570 c #SUrusnSV(ikl) = rusnSV(ikl) - zSlush ! [mm] OR [kg/m2]571 c #SURuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV)572 c #SUeta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn) !573 c #SU. *(1. - eta_SV(ikl,isn))) !574 c #SU. / max (ro_new , epsi ) !575 c #SUro__SV(ikl,isn) = ro_new !576 c #SUendif577 c #SUEND DO578 c #SUEND DO579 552 PorVol = 1. - ro__SV(ikl,isn) ! [--] 553 . *(1. - eta_SV(ikl,isn))/ ro_Ice ! 554 . - eta_SV(ikl,isn) ! 555 . *ro__SV(ikl,isn) / ro_Wat ! 556 PorVol = max(PorVol , zero ) ! 557 zWater = dzsnSV(ikl,isn) * PorVol * 1000. ! [mm] OR [kg/m2] 558 . * (1. -SWS_SV(ikl) ! 0 <=> freezing 559 . *(1 -min(1,iabs(isn-isnoSV(ikl))))) ! 1 <=> isn=isnoSV 560 zSlush = min(rusnSV(ikl) , zWater) ! [mm] OR [kg/m2] 561 ro_new =(dzsnSV(ikl,isn) * ro__SV(ikl,isn) ! 562 . +zSlush ) ! 563 . / max(dzsnSV(ikl,isn) , epsi ) ! 564 if(ro_new<ro_Ice+20) then ! MAX 940kg/m3 ! 565 rusnSV(ikl) = rusnSV(ikl) - zSlush ! [mm] OR [kg/m2] 566 RuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV) 567 eta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn) ! 568 . *(1. - eta_SV(ikl,isn))) ! 569 . / max (ro_new , epsi ) ! 570 ro__SV(ikl,isn) = ro_new ! 571 endif 572 END DO 573 END DO 574 END IF 580 575 581 576 C +--Impact of the Sublimation/Deposition on the Surface Mass Balance -
LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sop.F
r3792 r3900 1 2 3 1 subroutine SnOptP(jjtime) 4 2 … … 30 28 C | OUTPUT: albisv : Snow/Ice/Water/Soil Integrated Albedo [-] | 31 29 C | ^^^^^^ sEX_sv : Verticaly Integrated Extinction Coefficient | 30 C | DOPsnSV : Snow Optical diameter [m] | 32 31 C | | 33 32 C | Internal Variables: | … … 68 67 use VARySV 69 68 use VARtSV 70 USE surface_data, only: iflag_alb zenith69 USE surface_data, only: iflag_albcalc,correc_alb 71 70 72 71 IMPLICIT NONE … … 89 88 c #AG real agesno 90 89 91 integer isn ,ikl ,isn1 90 integer isn ,ikl ,isn1, i 92 91 real sbeta1,sbeta2,sbeta3,sbeta4,sbeta5 93 92 real AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK 94 93 real dalbed,dalbeS,dalbeW 95 real bsegal,czemax,csegal 94 real bsegal,czemax,csegal,csza 96 95 real RoFrez,DiffRo,SignRo,SnowOK,OpSqrt 97 96 real albSn1,albIc1,a_SnI1,a_SII1 … … 103 102 real albedo_old,albCor 104 103 real ro_ave,dz_ave,minalb 105 real thetazs,thetazs0,aap,bbp,ccp,alb0 ! parameters for zenoth angle effect at Dome C106 107 104 real l1min,l1max,l2min,l2max,l3min,l3max 105 real l6min(6), l6max(6), albSn6(6), a_SII6(6) 106 real lmintmp,lmaxtmp,albtmp 108 107 109 108 C +--Local DATA … … 138 137 data albmax /0.99 / ! Albedo max 139 138 140 C +-- For the computation of solar zenoth angle effect from Dome C data 141 C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 142 data aap/0.00016/,bbp/-0.017/,ccp/1.2/ 143 139 C +-- wavelength limits [m] for each broad band 140 141 data l1min/400.0e-9/,l1max/800.0e-9/ 142 data l2min/800.0e-9/,l2max/1500.0e-9/ 143 data l3min/1500.0e-9/,l3max/2800.0e-9/ 144 145 data l6min/185.0e-9,250.0e-9,400.0e-9, 146 . 690.0e-9,1190.0e-9,2380.0e-9/ 147 data l6max/250.0e-9,400.0e-9, 148 . 690.0e-9,1190.0e-9,2380.0e-9,4000.0e-9/ 144 149 145 150 … … 177 182 . *DFcdSV ))) 178 183 SnOpSV(ikl,isn) = max(zero,SnOpSV(ikl,isn)) 184 185 C + --For outputs (Etienne) 186 C + ------------------------ 187 DOPsnSV(ikl,isn)=SnOpSV(ikl,isn) 179 188 END DO 180 189 END DO 181 190 191 192 182 193 183 194 C +--Snow/Ice Albedo … … 191 202 DO ikl=1,knonv 192 203 193 204 isn = max(iun,isnoSV(ikl)) 194 205 195 206 SignRo = sign(unun, rocdSV - ro__SV(ikl,isn)) … … 200 211 cCA +--Correction of snow albedo for Antarctica/Greenland 201 212 cCA -------------------------------------------------- 202 albCor = 1. 213 214 215 albCor = correc_alb 203 216 c #GL albCor = 1.01 204 c #AC albCor = 1.01 205 217 c #AC albCor = 1.01 218 219 220 IF (iflag_albcalc .GE. 1) THEN ! Albedo calculation according to Kokhanovsky and Zege 2004 221 222 dalbed = 0.0 223 doptic=SnOpSV(ikl,isn) 224 csza=coszSV(ikl) 225 226 CALL albedo_kokhanovsky(l1min,l1max,csza,doptic,albSn1) 227 CALL albedo_kokhanovsky(l2min,l2max,csza,doptic,albSn2) 228 CALL albedo_kokhanovsky(l3min,l3max,csza,doptic,albSn3) 229 230 DO i=1,6 231 lmintmp=l6min(i) 232 lmaxtmp=l6max(i) 233 CALL albedo_kokhanovsky(lmintmp,lmaxtmp,csza,doptic,albtmp) 234 albSn6(i)=albtmp 235 ENDDO 236 237 238 ELSE ! Default calculation in SISVAT 239 240 ! Zenith Angle Correction (Segal et al., 1991, JAS 48, p.1025) 241 !--------------------------- (Wiscombe & Warren, dec1980, JAS , p.2723) 242 ! (Warren, 1982, RG , p. 81) 243 ! ------------------------------------------------- 244 245 dalbed = 0.0 246 247 csegal = max(czemax ,coszSV(ikl)) 248 c #cz dalbeS = ((bsegal+unun)/(unun+2.0*bsegal*csegal) 249 c #cz. - unun )*0.32 250 c #cz. / bsegal 251 c #cz dalbeS = max(dalbeS,zero) 252 c #cz dalbed = dalbeS * min(1,isnoSV(ikl)) 253 254 dalbeW =(0.64 - csegal )*0.0625 ! Warren 1982, RevGeo, fig.12b 255 ! 0.0625 = 5% * 1/0.8, p.81 256 ! 0.64 = cos(50) 257 dalbed = dalbeW * min(1,isnoSV(ikl)) 258 !------------------------------------------------------------------------- 259 206 260 albSn1 = 0.96-1.580*OpSqrt 207 261 albSn1 = max(albSn1,AlbMin) … … 219 273 albSn3 = min(albSn3*albCor,unun) 220 274 275 albSn6(1:3)=albSn1 276 albSn6(4:6)=albSn2 277 221 278 !snow albedo corection if wetsnow 222 279 c #GL albSn1 = albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn))) 223 280 c #GL albSn2 = albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn))) 224 281 c #GL albSn3 = albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn))) 282 283 ENDIF 284 225 285 226 286 albSno = So1dSV*albSn1 … … 237 297 albSn2 = SnowOK*albSn2+(1.0-SnowOK)*max(albSno,minalb) 238 298 albSn3 = SnowOK*albSn3+(1.0-SnowOK)*max(albSno,minalb) 299 albSn6(:) = SnowOK*albSn6(:)+(1.0-SnowOK)*max(albSno,minalb) 300 239 301 240 302 c + ro < roSdSV | min al > aI3dSV … … 304 366 a_SII3 = albWIc +(albSn3-albWIc) *SnownH 305 367 a_SII3 = min(a_SII3 ,albSn3) 306 368 369 DO i=1,6 370 a_SII6(i) = albWIc +(albSn6(i)-albWIc) *SnownH 371 a_SII6(i) = min(a_SII6(i) ,albSn6(i)) 372 ENDDO 373 307 374 cc #AG agesno = min(agsnSV(ikl,isn) ,AgeMax) 308 375 cc #AG a_SII1 = a_SII1 -0.175*agesno/AgeMax 309 376 C +... Impurities: Col de Porte Parameter. 310 377 311 312 ! Zenith Angle Correction (Segal et al., 1991, JAS 48, p.1025) 313 ! ----------------------- (Wiscombe & Warren, dec1980, JAS , p.2723) 314 ! (Warren, 1982, RG , p. 81) 315 ! (Vignon, PhD Thesis, pp 150, conditions at Dome C) 316 ! ------------------------------------------------- 317 318 319 dalbed = 0.0 320 321 if (iflag_albzenith .GE. 0) then 322 csegal = max(czemax ,coszSV(ikl)) 323 c #cz dalbeS = ((bsegal+unun)/(unun+2.0*bsegal*csegal) 324 c #cz. - unun )*0.32 325 c #cz. / bsegal 326 c #cz dalbeS = max(dalbeS,zero) 327 c #cz dalbed = dalbeS * min(1,isnoSV(ikl)) 328 329 dalbeW =(0.64 - csegal )*0.0625 ! Warren 1982, RevGeo, fig.12b 330 ! 0.0625 = 5% * 1/0.8, p.81 331 ! 0.64 = cos(50) 332 dalbed = dalbeW * min(1,isnoSV(ikl)) 333 334 335 ! Vignon PhD thesis, Dome C conditions 336 337 if (iflag_albzenith .EQ. 1) then 338 thetazs0=-bbp/(2.0*aap) 339 alb0=(bbp**2)/4./aap-(bbp**2)/2./aap+ccp 340 thetazs=max(acos(coszSV(ikl))/pi*180.0,thetazs0) ! in degrees 341 342 343 dalbeW = aap*(thetazs**2)+bbp*thetazs+ccp-alb0 344 dalbed = dalbeW * min(1,isnoSV(ikl)) 345 end if 346 347 348 end if 378 349 379 350 380 C +--Elsewhere Integrated Snow/Ice Albedo … … 368 398 alb3sv(ikl) = albssv(ikl) +(a_SII3-albssv(ikl))*SIcenH 369 399 alb3sv(ikl) = min(alb3sv(ikl) ,a_SII3) 370 400 371 401 albisv(ikl) = albssv(ikl) +(albSII-albssv(ikl))*SIcenH 372 402 albisv(ikl) = min(albisv(ikl) ,albSII) 373 403 404 DO i=1,6 405 alb6sv(ikl,i) = albssv(ikl) +(a_SII6(i)-albssv(ikl))*SIcenH 406 alb6sv(ikl,i) = min(alb6sv(ikl,i) ,a_SII6(i)) 407 ENDDO 408 374 409 375 410 C +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994 … … 382 417 alb3sv(ikl) = alb3sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH 383 418 . + dalbed * (1.-cld_SV(ikl)) 419 alb6sv(ikl,:) =alb6sv(ikl,:)+0.05 *(cld_SV(ikl)-0.5)*SIcenH 420 . + dalbed * (1.-cld_SV(ikl)) 384 421 albisv(ikl) = albisv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH 385 422 . + dalbed * (1.-cld_SV(ikl)) … … 397 434 alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 ! 33 % 398 435 . * (albedo_old-albisv(ikl)) / So3dSV 399 400 401 C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95% 436 alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 ! 16 % 437 . * (albedo_old-albisv(ikl)) / (So1dSV/3) 438 alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 ! 16 % 439 . * (albedo_old-albisv(ikl)) / (So2dSV/3) 440 441 442 C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95% 402 443 C + ----------------------------------------------------- 403 444 … … 410 451 alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 ! 33 % 411 452 . * (albedo_old-albisv(ikl)) / So3dSV 453 alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 ! 16 % 454 . * (albedo_old-albisv(ikl)) / (So1dSV/3) 455 alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 ! 16 % 456 . * (albedo_old-albisv(ikl)) / (So2dSV/3) 412 457 413 458 … … 436 481 alb2sv(ikl) = min(max(zero,alb2sv(ikl)),albmax) 437 482 alb3sv(ikl) = min(max(zero,alb3sv(ikl)),albmax) 438 483 484 DO i=1,6 485 alb6sv(ikl,i) = min(max(zero,alb6sv(ikl,i)),albmax) 486 ENDDO 439 487 END DO 440 488 … … 519 567 520 568 return 569 570 521 571 end 572 573 574 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 575 SUBROUTINE albedo_kokhanovsky(lambdamin,lambdamax, 576 . cossza,dopt,albint) 577 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 578 ! Authors: Hajar El Habchi El Fenniri, Etienne Vignon, Cecile Agosta 579 ! Ghislain Picard 580 ! Routine that calculates the integrated snow spectral albedo between 581 ! lambdamin and lambdamax following Kokhanisky and Zege 2004, 582 ! Scattering optics of snow, Applied Optics, Vol 43, No7 583 ! Code inspired from the snowoptics package of Ghislain Picard: 584 ! https://github.com/ghislainp/snowoptics 585 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 586 587 588 USE VARphy 589 590 IMPLICIT NONE 591 592 ! Inputs 593 !-------- 594 REAL lambdamin ! minimum wavelength for integration [m] 595 REAL lambdamax ! maximum wavelength for integration [m] 596 REAL cossza ! solar zenith angle cosinus 597 REAL dopt ! optical diameter [m] 598 599 !Outputs 600 !------- 601 REAL albint 602 603 ! Local Variables 604 !----------------- 605 606 REAL ropt,cosalb,norm,Pas 607 REAL SSA,alpha,gamm,R,cos30,alb30 608 INTEGER i 609 610 611 REAL B_amp ! amplification factor 612 PARAMETER (B_amp=1.6) 613 614 REAL g_asy ! asymetry factor 615 PARAMETER (g_asy=0.845) 616 617 INTEGER nlambda ! length of wavelength vector 618 PARAMETER(nlambda=200) 619 620 REAL lmin 621 PARAMETER(lmin=185.0E-9) 622 623 REAL lmax 624 PARAMETER(lmax=4000.0E-9) 625 626 REAL albmax 627 PARAMETER(albmax=0.99) 628 629 REAL wavelengths(nlambda) 630 REAL ni(nlambda) 631 632 DATA wavelengths / 1.85000000e-07, 2.04170854e-07, 633 . 2.23341709e-07, 2.42512563e-07, 634 . 2.61683417e-07, 2.80854271e-07, 3.00025126e-07, 3.19195980e-07, 635 . 3.38366834e-07, 3.57537688e-07, 3.76708543e-07, 3.95879397e-07, 636 . 4.15050251e-07, 4.34221106e-07, 4.53391960e-07, 4.72562814e-07, 637 . 4.91733668e-07, 5.10904523e-07, 5.30075377e-07, 5.49246231e-07, 638 . 5.68417085e-07, 5.87587940e-07, 6.06758794e-07, 6.25929648e-07, 639 . 6.45100503e-07, 6.64271357e-07, 6.83442211e-07, 7.02613065e-07, 640 . 7.21783920e-07, 7.40954774e-07, 7.60125628e-07, 7.79296482e-07, 641 . 7.98467337e-07, 8.17638191e-07, 8.36809045e-07, 8.55979899e-07, 642 . 8.75150754e-07, 8.94321608e-07, 9.13492462e-07, 9.32663317e-07, 643 . 9.51834171e-07, 9.71005025e-07, 9.90175879e-07, 1.00934673e-06, 644 . 1.02851759e-06, 1.04768844e-06, 1.06685930e-06, 1.08603015e-06, 645 . 1.10520101e-06, 1.12437186e-06, 1.14354271e-06, 1.16271357e-06, 646 . 1.18188442e-06, 1.20105528e-06, 1.22022613e-06, 1.23939698e-06, 647 . 1.25856784e-06, 1.27773869e-06, 1.29690955e-06, 1.31608040e-06, 648 . 1.33525126e-06, 1.35442211e-06, 1.37359296e-06, 1.39276382e-06, 649 . 1.41193467e-06, 1.43110553e-06, 1.45027638e-06, 1.46944724e-06, 650 . 1.48861809e-06, 1.50778894e-06, 1.52695980e-06, 1.54613065e-06, 651 . 1.56530151e-06, 1.58447236e-06, 1.60364322e-06, 1.62281407e-06, 652 . 1.64198492e-06, 1.66115578e-06, 1.68032663e-06, 1.69949749e-06, 653 . 1.71866834e-06, 1.73783920e-06, 1.75701005e-06, 1.77618090e-06, 654 . 1.79535176e-06, 1.81452261e-06, 1.83369347e-06, 1.85286432e-06, 655 . 1.87203518e-06, 1.89120603e-06, 1.91037688e-06, 1.92954774e-06, 656 . 1.94871859e-06, 1.96788945e-06, 1.98706030e-06, 2.00623116e-06, 657 . 2.02540201e-06, 2.04457286e-06, 2.06374372e-06, 2.08291457e-06, 658 . 2.10208543e-06, 2.12125628e-06, 2.14042714e-06, 2.15959799e-06, 659 . 2.17876884e-06, 2.19793970e-06, 2.21711055e-06, 2.23628141e-06, 660 . 2.25545226e-06, 2.27462312e-06, 2.29379397e-06, 2.31296482e-06, 661 . 2.33213568e-06, 2.35130653e-06, 2.37047739e-06, 2.38964824e-06, 662 . 2.40881910e-06, 2.42798995e-06, 2.44716080e-06, 2.46633166e-06, 663 . 2.48550251e-06, 2.50467337e-06, 2.52384422e-06, 2.54301508e-06, 664 . 2.56218593e-06, 2.58135678e-06, 2.60052764e-06, 2.61969849e-06, 665 . 2.63886935e-06, 2.65804020e-06, 2.67721106e-06, 2.69638191e-06, 666 . 2.71555276e-06, 2.73472362e-06, 2.75389447e-06, 2.77306533e-06, 667 . 2.79223618e-06, 2.81140704e-06, 2.83057789e-06, 2.84974874e-06, 668 . 2.86891960e-06, 2.88809045e-06, 2.90726131e-06, 2.92643216e-06, 669 . 2.94560302e-06, 2.96477387e-06, 2.98394472e-06, 3.00311558e-06, 670 . 3.02228643e-06, 3.04145729e-06, 3.06062814e-06, 3.07979899e-06, 671 . 3.09896985e-06, 3.11814070e-06, 3.13731156e-06, 3.15648241e-06, 672 . 3.17565327e-06, 3.19482412e-06, 3.21399497e-06, 3.23316583e-06, 673 . 3.25233668e-06, 3.27150754e-06, 3.29067839e-06, 3.30984925e-06, 674 . 3.32902010e-06, 3.34819095e-06, 3.36736181e-06, 3.38653266e-06, 675 . 3.40570352e-06, 3.42487437e-06, 3.44404523e-06, 3.46321608e-06, 676 . 3.48238693e-06, 3.50155779e-06, 3.52072864e-06, 3.53989950e-06, 677 . 3.55907035e-06, 3.57824121e-06, 3.59741206e-06, 3.61658291e-06, 678 . 3.63575377e-06, 3.65492462e-06, 3.67409548e-06, 3.69326633e-06, 679 . 3.71243719e-06, 3.73160804e-06, 3.75077889e-06, 3.76994975e-06, 680 . 3.78912060e-06, 3.80829146e-06, 3.82746231e-06, 3.84663317e-06, 681 . 3.86580402e-06, 3.88497487e-06, 3.90414573e-06, 3.92331658e-06, 682 . 3.94248744e-06, 3.96165829e-06, 3.98082915e-06, 4.00000000e-06/ 683 684 685 DATA ni /7.74508407e-10, 7.74508407e-10, 686 . 7.74508407e-10, 7.74508407e-10, 687 . 7.74508407e-10, 7.74508407e-10, 7.74508407e-10, 7.74508407e-10, 688 . 6.98381122e-10, 6.23170274e-10, 5.97655992e-10, 5.84106004e-10, 689 . 5.44327597e-10, 5.71923510e-10, 6.59723827e-10, 8.05183870e-10, 690 . 1.03110161e-09, 1.36680386e-09, 1.85161253e-09, 2.56487751e-09, 691 . 3.56462855e-09, 4.89450926e-09, 6.49252022e-09, 9.62029335e-09, 692 . 1.32335015e-08, 1.75502184e-08, 2.19240625e-08, 3.03304156e-08, 693 . 4.07715972e-08, 5.00414911e-08, 7.09722331e-08, 1.00773751e-07, 694 . 1.31427409e-07, 1.42289041e-07, 1.49066787e-07, 2.01558515e-07, 695 . 2.99106105e-07, 4.03902086e-07, 4.54292169e-07, 5.21906983e-07, 696 . 6.27643362e-07, 9.43955678e-07, 1.33464494e-06, 1.97278315e-06, 697 . 2.31801329e-06, 2.20584676e-06, 1.85568138e-06, 1.73395193e-06, 698 . 1.73101406e-06, 1.91333936e-06, 2.26413019e-06, 3.23959718e-06, 699 . 4.94316963e-06, 6.89378896e-06, 1.02237444e-05, 1.21439656e-05, 700 . 1.31567585e-05, 1.33448288e-05, 1.32000000e-05, 1.31608040e-05, 701 . 1.33048369e-05, 1.40321464e-05, 1.51526244e-05, 1.80342858e-05, 702 . 3.82875736e-05, 1.07325259e-04, 2.11961637e-04, 3.82008054e-04, 703 . 5.30897470e-04, 5.29244735e-04, 4.90876605e-04, 4.33905427e-04, 704 . 3.77795349e-04, 3.17633099e-04, 2.81078564e-04, 2.57579485e-04, 705 . 2.42203100e-04, 2.23789060e-04, 2.04306870e-04, 1.87909255e-04, 706 . 1.73117146e-04, 1.61533186e-04, 1.53420328e-04, 1.47578033e-04, 707 . 1.42334776e-04, 1.35691466e-04, 1.30495414e-04, 1.36065123e-04, 708 . 1.70928821e-04, 2.66389730e-04, 4.80957955e-04, 8.25041961e-04, 709 . 1.21654792e-03, 1.50232875e-03, 1.62316078e-03, 1.61649750e-03, 710 . 1.53736801e-03, 1.42343711e-03, 1.24459117e-03, 1.02388611e-03, 711 . 7.89112523e-04, 5.97204264e-04, 4.57152413e-04, 3.62341259e-04, 712 . 2.99128332e-04, 2.57035569e-04, 2.26992203e-04, 2.07110355e-04, 713 . 2.05835688e-04, 2.25108810e-04, 2.64262893e-04, 3.23594011e-04, 714 . 3.93061117e-04, 4.62789970e-04, 5.19664416e-04, 5.59739628e-04, 715 . 5.93476084e-04, 6.22797885e-04, 6.57484833e-04, 6.92849600e-04, 716 . 7.26584901e-04, 7.56604648e-04, 7.68009488e-04, 7.65579073e-04, 717 . 7.50526164e-04, 7.39809972e-04, 7.55622847e-04, 8.05099514e-04, 718 . 9.67279246e-04, 1.16281559e-03, 1.42570247e-03, 2.04986585e-03, 719 . 2.93971170e-03, 4.49827711e-03, 7.32537532e-03, 1.18889734e-02, 720 . 1.85851805e-02, 2.86242532e-02, 4.34131035e-02, 6.37828307e-02, 721 . 9.24145850e-02, 1.35547945e-01, 1.94143245e-01, 2.54542814e-01, 722 . 3.02282024e-01, 3.42214181e-01, 3.85475065e-01, 4.38000000e-01, 723 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 724 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 725 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 726 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 727 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 728 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 729 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 730 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 731 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 732 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 733 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 734 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 735 . 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01/ 736 737 738 Pas =(lmax-lmin)/nlambda 739 ropt = dopt/2.0 740 SSA = 3.0/(rhoIce*ropt) 741 cos30 = cos(30.0/180.0*pi) 742 743 744 albint=0. 745 norm=0. 746 747 DO i = 1,nlambda 748 gamm = 4.0 * pi * ni(i) / wavelengths(i) 749 cosalb = 2.0 / (SSA * rhoice) * B_amp * gamm 750 alpha = 16. / 3 * cosalb / (1.0 - g_asy) 751 R = exp(-(alpha**0.5) * 3.0 / 7.0 * (1.0 + 2.0 * cossza)) 752 alb30 = exp(-(alpha**0.5)* 3.0 / 7.0 * (1.0 + 20 * cos30)) 753 754 IF ((wavelengths(i).GE.lambdamin).AND. 755 . (wavelengths(i).LT.lambdamax)) THEN 756 albint=albint+R*Pas ! rectangle integration -> can be improved with trapezintegration 757 norm=norm+Pas 758 ENDIF 759 760 END DO 761 762 albint=max(0.,min(albint/max(norm,1E-10),albmax)) 763 764 END 765 -
LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_tso.F
r3792 r3900 132 132 133 133 integer nt_srf,it_srf,itEuBk ! HL: Surface Scheme 134 parameter(nt_srf=10) !134 parameter(nt_srf=10) ! 10 before 135 135 real agpsrf,xgpsrf,dt_srf,dt_ver ! 136 136 real etaBAK(knonv) ! … … 153 153 C + ! including Snow Melt Energy 154 154 155 155 C +-- Initilialisation of local arrays 156 C + ================================ 157 DO ikl=1,knonv 158 159 mu_sno(ikl)=0. 160 mu__dz(ikl,:)=0. 161 dtC_sv(ikl,:)=0. 162 IRs__D(ikl)=0. 163 dIRsdT(ikl)=0. 164 f_HSHL(ikl)=0. 165 dRidTs(ikl)=0. 166 HS___D(ikl)=0. 167 f___HL(ikl)=0. 168 HL___D(ikl)=0. 169 TSurf0(ikl)=0. 170 qsatsg(ikl)=0. 171 dqs_dT(ikl)=0. 172 Psi(ikl)=0. 173 RHuSol(ikl)=0. 174 Diag_A(ikl,:)=0. 175 Diag_B(ikl,:)=0. 176 Diag_C(ikl,:)=0. 177 Term_D(ikl,:)=0. 178 Aux__P(ikl,:)=0. 179 Aux__Q(ikl,:)=0. 180 etaBAK(ikl)=0. 181 etaNEW(ikl)=0. 182 etEuBk(ikl)=0. 183 fac_dt(ikl)=0. 184 faceta(ikl)=0. 185 PsiArg(ikl)=0. 186 SHuSol(ikl)=0. 187 188 END DO 189 190 156 191 157 192 C +--Heat Conduction Coefficient (zero in the Layers over the highest one) … … 336 371 C +--Snow highest Layer (dummy!) 337 372 C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^ 338 isl= min(isnoSV(1)+1,nsno) 339 DO ikl=1,knonv 373 374 !EV!isl= min(isnoSV(1)+1,nsno) 375 376 DO ikl=1,knonv 377 ! EV try to calculate isl at the ikl grid point 378 isl= min(isnoSV(ikl)+1,nsno) 379 340 380 Elem_A = dtC_sv(ikl,isl) *mu__dz(ikl,isl) 341 381 Elem_C = 0. … … 384 424 c . / den_qs ! 385 425 c qsatsg(ikl) = .0038 * exp(arg_qs) ! 386 387 426 ! sp = (pst_SV(ikl) + ptopSV) * 10. 388 427 389 sp=ps__SV(ikl) 428 !sp=ps__SV(ikl) 429 ! Etienne: in the formula herebelow sp should be in hPa, not 430 ! in Pa so I divide by 100. 431 sp=ps__SV(ikl)/100. 390 432 psat_ice = 6.1070 * exp(6150. *(1./273.16 - 391 433 . 1./TsisSV(ikl,isl))) … … 399 441 qsatsg(ikl) = 0.622 * psat_wat / (sp - 0.378 * psat_wat) 400 442 endif 443 QsT_SV(ikl)=qsatsg(ikl) 401 444 402 445 c dqs_dT(ikl) = qsatsg(ikl)* 4099.2 /(den_qs *den_qs)! 403 446 fac_dt(ikl) = f_HSHL(ikl)/(ro_Wat * dz_dSV(0)) ! 404 447 END DO 448 449 405 450 406 451 C +--Surface: Latent Heat Flux: Surface Relative Humidity … … 410 455 . /( 1.0-xgpsrf**nt_srf) ! 411 456 dt_srf = agpsrf ! 412 dt_ver = 0. ! 457 dt_ver = 0. 458 413 459 DO ikl=1,knonv 414 isl = isnoSV(ikl) ! 460 isl = isnoSV(ikl) 461 ist = max(0,isotSV(ikl)-100*isnoSV(ikl))! 0 if H2O 462 ist__s = min(1,ist) 415 463 etaBAK(ikl) = max(epsi,eta_SV(ikl ,isl)) ! 416 464 etaNEW(ikl) = etaBAK(ikl) ! 417 465 etEuBk(ikl) = etaNEW(ikl) ! 418 END DO ! 466 END DO 467 468 if(ist__s==1) then ! to reduce computer time 469 ! 419 470 DO it_srf=1,nt_srf ! 420 471 dt_ver = dt_ver +dt_srf ! … … 458 509 END DO ! 459 510 dt_srf = dt_srf * xgpsrf ! 460 END DO ! 511 END DO 512 513 514 endif ! 461 515 462 516 C +--Surface: Latent Heat Flux: Soil/Water Surface Contributions … … 579 633 580 634 END DO 635 636 581 637 582 638 C +--Temperature Limits (avoids problems in case of no Snow Layers) … … 584 640 DO ikl= 1,knonv 585 641 isl = isnoSV(ikl) 586 dTSurf = TsisSV(ikl,isl) - TSurf0(ikl) 642 643 dTSurf = TsisSV(ikl,isl) - TSurf0(ikl) 587 644 TsisSV(ikl,isl) = TSurf0(ikl) + sign(1.,dTSurf) ! 180.0 dgC/hr 588 645 . * min(abs(dTSurf),5.e-2*dt__SV) ! =0.05 dgC/s … … 602 659 C +--Update Surface Fluxes 603 660 C + ======================== 604 661 662 663 605 664 DO ikl= 1,knonv 606 665 isl = isnoSV(ikl) … … 613 672 END DO 614 673 615 616 617 674 return 618 675 end -
LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_zsn.F
r3792 r3900 52 52 use VARxSV 53 53 use VARySV 54 use surface_data, only: ok_zsn_ii 54 55 55 56 IMPLICIT NONE … … 716 717 END DO 717 718 719 720 C +--Search new Ice/Snow Interface (option II in MAR) 721 C + =============================================== 722 723 IF (ok_zsn_ii) THEN 724 725 DO ikl=1,knonv 726 iiceSV(ikl) = 0 727 END DO 728 729 DO ikl=1,knonv 730 DO isn=1,isnoSV(ikl) 731 OK_ICE = max(zero,sign(unun,ro__SV(ikl,isn)-ro_ice+20.)) 732 . * max(zero,sign(unun,dzsnSV(ikl,isn)-epsi)) 733 iiceSV(ikl) = (1.-OK_ICE) *iiceSV(ikl) 734 . + OK_ICE *isn 735 END DO 736 END DO 737 738 END IF 718 739 719 740 return -
LMDZ6/trunk/libf/phylmd/inlandsis/surf_inlandsis_mod.F90
r3792 r3900 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 -
LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
r3888 r3900 182 182 debut, lafin, & 183 183 rlon, rlat, rugoro, rmu0, & 184 zsig, lwdown_m, pphi,cldt, &184 lwdown_m, cldt, & 185 185 rain_f, snow_f, solsw_m, solswfdiff_m, sollw_m, & 186 186 gustiness, & … … 277 277 ! z0m, z0h ----input-R- longeur de rugosite (en m) 278 278 ! Martin 279 ! zsig-----input-R- slope280 279 ! cldt-----input-R- total cloud fraction 281 ! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol)282 280 ! Martin 283 281 ! … … 315 313 USE print_control_mod, ONLY : prt_level,lunout 316 314 USE ioipsl_getin_p_mod, ONLY : getin_p 317 use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal 315 use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, zsig, zmea 318 316 use phys_output_var_mod, only: dter, dser, tkt, tks, taur, sss 319 317 #ifdef CPP_XIOS … … 322 320 use netcdf, only: missing_val => nf90_fill_real 323 321 #endif 322 323 324 324 325 325 326 IMPLICIT NONE … … 359 360 REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! sub-surface fraction 360 361 ! Martin 361 REAL, DIMENSION(klon), INTENT(IN) :: zsig ! slope362 362 REAL, DIMENSION(klon), INTENT(IN) :: lwdown_m ! downward longwave radiation at mean s 363 363 REAL, DIMENSION(klon), INTENT(IN) :: gustiness ! gustiness 364 364 365 365 REAL, DIMENSION(klon), INTENT(IN) :: cldt ! total cloud fraction 366 REAL, DIMENSION(klon,klev), INTENT(IN) :: pphi ! geopotential (m2/s2)367 ! Martin368 366 369 367 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012 … … 873 871 REAL, DIMENSION(klon) :: ytoice 874 872 REAL, DIMENSION(klon) :: ysnowhgt, yqsnow, ysissnow, yrunoff 873 REAL, DIMENSION(klon) :: yzmea 875 874 REAL, DIMENSION(klon) :: yzsig 876 REAL, DIMENSION(klon,klev) :: ypphi877 875 REAL, DIMENSION(klon) :: ycldt 878 876 REAL, DIMENSION(klon) :: yrmu0 … … 1054 1052 ysnowhgt = 0.0; yqsnow = 0.0 ; yrunoff = 0.0 ; ytoice =0.0 1055 1053 yalb3_new = 0.0 ; ysissnow = 0.0 1056 y pphi = 0.0 ; ycldt = 0.0 ; yrmu0 = 0.01054 ycldt = 0.0 ; yrmu0 = 0.0 1057 1055 ! Martin 1058 1056 … … 1394 1392 ywindsp(j) = windsp(i,nsrf) 1395 1393 !>jyg 1396 ! Martin 1394 ! Martin and Etienne 1395 yzmea(j) = zmea(i) 1397 1396 yzsig(j) = zsig(i) 1398 1397 ycldt(j) = cldt(i) … … 1529 1528 ! 1530 1529 !**************************************************************************************** 1530 1531 1531 1532 1532 !!! jyg le 07/02/2012 … … 1579 1579 speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2) 1580 1580 ENDDO 1581 CALL cdrag(knon, nsrf, & 1581 1582 1583 CALL cdrag(knon, nsrf, & 1582 1584 speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),& 1583 1585 yts_x, yqsurf_x, yz0m, yz0h, & … … 2040 2042 CALL surf_landice(itap, dtime, knon, ni, & 2041 2043 rlon, rlat, debut, lafin, & 2042 yrmu0, ylwdown, yalb, ypphi(:,1), &2044 yrmu0, ylwdown, yalb, zgeo1, & 2043 2045 ysolsw, ysollw, yts, ypplay(:,1), & 2044 2046 !!jyg ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& … … 2050 2052 ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, & 2051 2053 ytsurf_new, y_dflux_t, y_dflux_q, & 2052 yz sig, ycldt, &2054 yzmea, yzsig, ycldt, & 2053 2055 ysnowhgt, yqsnow, ytoice, ysissnow, & 2054 2056 yalb3_new, yrunoff, & … … 2840 2842 sss(ni(:knon)) = ysss(:knon) 2841 2843 end if 2844 2842 2845 2843 2846 !**************************************************************************************** … … 3253 3256 ENDDO 3254 3257 !!! 3255 3258 3256 3259 ! 3257 3260 ! Incrementer la temperature du sol -
LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
r3893 r3900 1054 1054 CALL histwrite_phy(o_tauy, zx_tmp_fi2d) 1055 1055 1056 IF (landice_opt .GE. 1) THEN 1057 CALL histwrite_phy(o_snowsrf, snow_o) 1058 CALL histwrite_phy(o_qsnow, qsnow) 1059 CALL histwrite_phy(o_snowhgt,snowhgt) 1060 CALL histwrite_phy(o_toice,to_ice) 1061 CALL histwrite_phy(o_sissnow,sissnow) 1062 CALL histwrite_phy(o_runoff,runoff) 1063 CALL histwrite_phy(o_albslw3,albsol3_lic) 1064 ENDIF 1056 ! Etienne: test sorties pour compil sur JZ 1057 ! IF (landice_opt .GE. 1) THEN 1058 ! CALL histwrite_phy(o_snowsrf, snow_o) 1059 ! CALL histwrite_phy(o_qsnow, qsnow) 1060 ! CALL histwrite_phy(o_snowhgt,snowhgt) 1061 ! CALL histwrite_phy(o_toice,to_ice) 1062 ! CALL histwrite_phy(o_sissnow,sissnow) 1063 ! CALL histwrite_phy(o_runoff,runoff) 1064 ! CALL histwrite_phy(o_albslw3,albsol3_lic) 1065 ! ENDIF 1065 1066 1066 1067 DO nsrf = 1, nbsrf -
LMDZ6/trunk/libf/phylmd/physiq_mod.F90
r3888 r3900 2537 2537 debut, lafin, & 2538 2538 longitude_deg, latitude_deg, rugoro, zrmu0, & 2539 zsig, sollwdown, pphi, cldt, &2539 sollwdown, cldt, & 2540 2540 rain_fall, snow_fall, solsw, solswfdiff, sollw, & 2541 2541 gustiness, & -
LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90
r3792 r3900 19 19 tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, & 20 20 tsurf_new, dflux_s, dflux_l, & 21 slope, cloudf, &21 alt, slope, cloudf, & 22 22 snowhgt, qsnow, to_ice, sissnow, & 23 23 alb3, runoff, & … … 25 25 26 26 USE dimphy 27 USE surface_data, ONLY : type_ocean, calice, calsno, landice_opt, n_dtis28 USE fonte_neige_mod, ONLY : fonte_neige, run_off_lic27 USE surface_data, ONLY : type_ocean, calice, calsno, landice_opt, iflag_albcalc 28 USE fonte_neige_mod, ONLY : fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global 29 29 USE cpl_mod, ONLY : cpl_send_landice_fields 30 30 USE calcul_fluxs_mod … … 75 75 REAL, DIMENSION(klon), INTENT(IN) :: albedo !mean albedo 76 76 REAL, DIMENSION(klon), INTENT(IN) :: pphi1 77 REAL, DIMENSION(klon), INTENT(IN) :: alt !mean altitude of the grid box 77 78 REAL, DIMENSION(klon), INTENT(IN) :: slope !mean slope in grid box 78 79 REAL, DIMENSION(klon), INTENT(IN) :: cloudf !total cloud fraction … … 115 116 REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay, ustar 116 117 INTEGER :: i,j,nt 117 118 REAL, DIMENSION(klon) :: fqfonte,ffonte 118 119 REAL, DIMENSION(klon) :: emis_new !Emissivity 119 120 REAL, DIMENSION(klon) :: swdown,lwdown 120 REAL, DIMENSION(klon) :: precip_snow_adv, snow_adv !Snow Drift precip./advection 121 REAL, DIMENSION(klon) :: zsl_height, wind_velo !surface layer height, wind spd 121 REAL, DIMENSION(klon) :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis) 122 REAL, DIMENSION(klon) :: erod !erosion of surface snow (flux, kg/m2/s like evap) 123 REAL, DIMENSION(klon) :: zsl_height, wind_velo !surface layer height, wind spd 122 124 REAL, DIMENSION(klon) :: dens_air, snow_cont_air !air density; snow content air 123 125 REAL, DIMENSION(klon) :: alb_soil !albedo of underlying ice … … 132 134 133 135 134 !albedo SB >>> 135 real,dimension(klon) :: alb1,alb2 136 !albedo SB <<< 137 136 REAL,DIMENSION(klon) :: alb1,alb2 137 REAL, DIMENSION (klon,6) :: alb6 138 138 ! End definition 139 139 !**************************************************************************************** … … 230 230 AcoefH, AcoefQ, BcoefH, BcoefQ, cdragh, & 231 231 run_off_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, & 232 tsurf_new, alb1, alb2, alb3, & 233 emis_new, z0m, qsurf) 232 tsurf_new, alb1, alb2, alb3, emis_new, z0m, qsurf) 234 233 z0h(1:knon)=z0m(1:knon) ! en attendant mieux 235 234 … … 278 277 swdown(:) = 0.0 279 278 lwdown(:) = 0.0 280 snow_adv(:) = 0. ! no snow blown in for now 281 snow_cont_air(:) = 0. 279 snow_cont_air(:) = 0. ! the snow content in air is not a prognostic variable of the model 282 280 alb_soil(:) = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set 283 281 ustar(:) = 0. … … 296 294 297 295 298 ! Subtimestepping 299 300 dtis=dtime/n_dtis 301 302 DO nt=1,n_dtis 303 304 IF (lafin .and. nt.eq.n_dtis) THEN 296 297 dtis=dtime 298 299 IF (lafin) THEN 305 300 lafin_is=.true. 306 301 END IF 307 302 308 !PRINT*,'RENTRE DANS INLANDSIS','itime',itime,'dtime',dtime,'dtis',dtis 309 CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is, & 310 rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, & 311 precip_rain, precip_snow, precip_snow_adv, snow_adv, & 312 zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, & 313 rugoro, snow_cont_air, alb_soil, slope, cloudf, & 314 radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice,sissnow, agesno, & 303 CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,& 304 rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow, & 305 zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,& 306 rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, & 307 radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno, & 315 308 AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, & 316 run_off_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, & 317 tsurf_new, alb1, alb2, alb3, & 318 emis_new, z0m, z0h, qsurf) 319 320 debut_is=.false. 321 322 END DO 309 run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, & 310 tsurf_new, alb1, alb2, alb3, alb6, & 311 emis_new, z0m, z0h, qsurf) 312 313 debut_is=.false. 314 315 316 ! Treatment of snow melting and calving 317 318 ! for consistency with standard LMDZ, add calving to run_off_lic 319 run_off_lic(:)=run_off_lic(:) + to_ice(:) 320 321 DO i = 1, knon 322 ffonte_global(knindex(i),is_lic) = ffonte(i) 323 fqfonte_global(knindex(i),is_lic) = fqfonte(i)! net melting= melting - refreezing 324 fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux 325 runofflic_global(knindex(i)) = run_off_lic(i) 326 ENDDO 327 ! Here, we assume that the calving term is equal to the to_ice term 328 ! (no ice accumulation) 323 329 324 330 … … 419 425 420 426 421 422 423 424 427 425 428 END IF ! landice_opt … … 476 479 alb_dir(1:knon,5)=alb2(1:knon) 477 480 alb_dir(1:knon,6)=alb2(1:knon) 481 482 IF ((landice_opt .EQ.2) .AND. (iflag_albcalc .EQ. 2)) THEN 483 alb_dir(1:knon,1)=alb6(1:knon,1) 484 alb_dir(1:knon,2)=alb6(1:knon,2) 485 alb_dir(1:knon,3)=alb6(1:knon,3) 486 alb_dir(1:knon,4)=alb6(1:knon,4) 487 alb_dir(1:knon,5)=alb6(1:knon,5) 488 alb_dir(1:knon,6)=alb6(1:knon,6) 489 ENDIF 490 478 491 end select 479 492 alb_dif=alb_dir 480 493 !albedo SB <<< 481 494 482 483 495 496 484 497 485 498 END SUBROUTINE surf_landice -
LMDZ6/trunk/libf/phylmd/surface_data.F90
r3792 r3900 29 29 ! FOR INLANDSIS: 30 30 !=============== 31 32 INTEGER, SAVE :: landice_opt ! 1 for coupling with SISVAT, 2 for coupling with INLANDSIS 31 32 ! 1 for coupling with SISVAT, 2 for coupling with INLANDSIS and number of subtimesteps for INLANDSIS 33 INTEGER, SAVE :: landice_opt ! 1 for coupling with SISVAT, 2 for coupling with INLANDSIS 33 34 !$OMP THREADPRIVATE(landice_opt) 34 35 35 INTEGER, SAVE :: iflag_tsurf_inlandsis ! 0 SISVAT method, 1 LMDZ method 36 !$OMP THREADPRIVATE(iflag_tsurf_inlandsis) 36 ! temperature calculation options within the soil and at the surface 37 INTEGER, SAVE :: iflag_tsurf_inlandsis,iflag_temp_inlandsis 38 !$OMP THREADPRIVATE(iflag_tsurf_inlandsis,iflag_temp_inlandsis) 37 39 38 INTEGER, SAVE :: iflag_albzenith ! dependency of albedo to zenith angle 39 !$OMP THREADPRIVATE(iflag_albzenith) 40 41 INTEGER, SAVE :: n_dtis ! number of subtimesteps for INLANDSIS 42 !$OMP THREADPRIVATE(n_dtis) 40 ! flags for albedo and roughness calc. 41 INTEGER, SAVE :: iflag_albcalc,iflag_z0m_snow 42 !$OMP THREADPRIVATE(iflag_albcalc,iflag_z0m_snow) 43 43 44 44 ! with or without snow module/ blowing snow, ascii outfile 45 LOGICAL, SAVE 45 LOGICAL, SAVE :: SnoMod,BloMod,ok_outfor 46 46 !$OMP THREADPRIVATE(SnoMod,BloMod,ok_outfor) 47 47 48 ! activate slush, korlyakov snow density, RN z0h calc. 49 LOGICAL, SAVE :: is_ok_slush,is_ok_density_kotlyakov,is_ok_z0h_rn 50 !$OMP THREADPRIVATE(is_ok_slush,is_ok_density_kotlyakov,is_ok_z0h_rn) 51 52 ! activate detection snow/ice layers and option XF discrtet/option runoff AC 53 LOGICAL, SAVE :: ok_zsn_ii,discret_xf,opt_runoff_ac 54 !$OMP THREADPRIVATE(ok_zsn_ii,discret_xf, opt_runoff_ac) 55 56 ! value of z0m snow when prescribed and albedo correction term 57 REAL, SAVE :: prescribed_z0m_snow,correc_alb 58 !$OMP THREADPRIVATE(prescribed_z0m_snow, correc_alb) 59 60 ! value of sphericity [0-99] and snow grain size [e-4m] for polar buffer snow 61 ! layer 62 REAL, SAVE :: buf_sph_pol,buf_siz_pol 63 !$OMP THREADPRIVATE(buf_sph_pol,buf_siz_pol) 64 65 66 48 67 END MODULE surface_data
Note: See TracChangeset
for help on using the changeset viewer.