- Timestamp:
- Jul 28, 2025, 7:23:15 PM (6 days ago)
- Location:
- LMDZ6/branches/contrails
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/ocean_forced_mod.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/contrails
- Property svn:mergeinfo changed
/LMDZ6/trunk merged: 5654-5683,5685-5690,5692-5715,5718-5721,5726-5727,5729,5744-5761,5763-5778,5780,5785-5789
- Property svn:mergeinfo changed
-
LMDZ6/branches/contrails/libf/phylmd/ocean_forced_mod.F90
r5618 r5791 22 22 radsol, snow, agesno, & 23 23 qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & 24 tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa & 24 tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa, & 25 dthetadz300,pctsrf,Ampl & 25 26 #ifdef ISO 26 27 ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, & … … 39 40 USE mod_grid_phy_lmdz 40 41 USE indice_sol_mod 42 USE surface_data, ONLY : iflag_leads 41 43 USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o 42 44 use config_ocean_skin_m, only: activate_ocean_skin … … 69 71 REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in 70 72 real, intent(in):: rhoa(:) ! (knon) density of moist air (kg / m3) 73 !GG 74 REAL, DIMENSION(klon), INTENT(IN) :: dthetadz300 75 REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf 76 ! 71 77 72 78 #ifdef ISO … … 93 99 REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new 94 100 REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l 95 REAL, intent(out):: sens_prec_liq(:) ! (knon) 101 REAL, INTENT(out):: sens_prec_liq(:) ! (knon) 102 !GG 103 REAL, DIMENSION(klon), INTENT(OUT) :: Ampl 104 ! 96 105 97 106 #ifdef ISO … … 110 119 REAL, DIMENSION(knon) :: sens_prec_sol 111 120 REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol 121 ! GG 122 REAL, DIMENSION(klon) :: l_CBL, sicfra 123 ! 112 124 #ifdef ISO 113 125 REAL, PARAMETER :: t_coup = 273.15 … … 204 216 enddo 205 217 218 !GG 219 if (iflag_leads == 1) then 220 l_CBL = -52381.*dthetadz300 + 3008.1 221 Ampl = 6.012e-08*l_CBL**2 - 4.036e-04*l_CBL + 1.4979 222 WHERE(Ampl(:)>1.2) Ampl(:)=1.2 223 sicfra(:)=pctsrf(:,is_sic)/(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 224 WHERE(pctsrf(:,is_sic)+pctsrf(:,is_oce)<EPSFRA) sicfra(:)=0. 225 WHERE(sicfra<0.7) Ampl(:)=1. 226 WHERE((sicfra>0.7).and.(sicfra<0.9)) Ampl=((sicfra-0.7)/0.2)*Ampl+((0.9-sicfra)/0.2) 227 fluxsens=Ampl*fluxsens 228 dflux_s=Ampl*dflux_s 229 endif 230 206 231 207 232 ! - Flux calculation at first modele level for U and V … … 244 269 AcoefH, AcoefQ, BcoefH, BcoefQ, & 245 270 AcoefU, AcoefV, BcoefU, BcoefV, & 246 ps, u1, v1, gustiness, & 271 !GG ps, u1, v1, gustiness, & 272 ps, u1, v1, gustiness, pctsrf, & 273 !GG 247 274 radsol, snow, qsol, agesno, tsoil, & 248 275 qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & 249 tsurf_new, dflux_s, dflux_l, rhoa & 276 !GG tsurf_new, dflux_s, dflux_l, rhoa) 277 tsurf_new, dflux_s, dflux_l, rhoa, swnet, hice, tice, bilg_cumul, & 278 fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, & 279 dtice_melt, dtice_snow2sic & 280 !GG 250 281 #ifdef ISO 251 282 ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, & … … 261 292 USE geometry_mod, ONLY: longitude,latitude 262 293 USE calcul_fluxs_mod 263 USE surface_data, ONLY : calice, calsno 294 !GG USE surface_data, ONLY : calice, calsno 295 USE surface_data, ONLY : calice, calsno, iflag_seaice, iflag_seaice_alb, & 296 sice_cond, sisno_cond, sisno_den, sisno_min, sithick_min, sisno_wfact, & 297 amax_s,amax_n, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, & 298 si_pen_frac, si_pen_ext, fseaN, fseaS, iflag_leads 299 300 USE geometry_mod, ONLY: longitude,latitude,latitude_deg 301 !GG 264 302 USE limit_read_mod 265 303 USE fonte_neige_mod, ONLY : fonte_neige … … 298 336 REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness 299 337 real, intent(in):: rhoa(:) ! (knon) density of moist air (kg / m3) 338 !GG 339 REAL, DIMENSION(klon), INTENT(IN) :: swnet 340 REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf 341 !GG 300 342 #ifdef ISO 301 343 REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow … … 311 353 REAL, DIMENSION(klon), INTENT(INOUT) :: agesno 312 354 REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil 355 !GG 356 REAL, DIMENSION(klon), INTENT(INOUT) :: hice 357 REAL, DIMENSION(klon), INTENT(INOUT) :: tice 358 REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul 359 REAL, DIMENSION(klon), INTENT(INOUT) :: fcds 360 REAL, DIMENSION(klon), INTENT(INOUT) :: fcdi 361 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_growth 362 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_melt 363 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt 364 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_snow2sic 365 REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_melt 366 REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_snow2sic 367 !GG 313 368 #ifdef ISO 314 369 REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsnow … … 342 397 REAL, DIMENSION(knon) :: sens_prec_liq, sens_prec_sol 343 398 REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol 399 !GG 400 INTEGER :: ki 401 INTEGER :: cpl_pas 402 REAL, DIMENSION(klon) :: bilg, fsic, f_bot 403 REAL, PARAMETER :: latent_ice = 334.0e3 404 REAL, PARAMETER :: rau_ice = 917.0 405 REAL, PARAMETER :: kice=2.2 406 REAL :: f_cond, f_swpen, f_cond_s, f_cond_i 407 REAL :: ustar, uscap, ustau 408 ! for snow/ice albedo: 409 REAL :: alb_snow, alb_ice, alb_pond 410 REAL :: frac_snow, frac_ice, frac_pond 411 REAL :: z1_i, z2_i, z1_s, zlog ! height parameters 412 ! for ice melt / freeze 413 REAL :: e_melt, snow_evap, h_test 414 ! dhsic, dfsic change in ice mass, fraction. 415 REAL :: dhsic, dfsic, frac_mf 416 REAL :: fsea, amax 417 REAL :: hice_i, tice_i, fsic_new 418 ! snow and ice physical characteristics: 419 REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp 420 REAL, PARAMETER :: t_melt=273.15 ! melting ice temp 421 REAL :: sno_den!=sisno_den !mean snow density, kg/m3 422 REAL, PARAMETER :: ice_den=917. ! ice density 423 REAL, PARAMETER :: sea_den=1025. ! sea water density 424 REAL :: ice_cond!=sice_cond*ice_den !conductivity of ice 425 REAL :: sno_cond!=sisno_cond*sno_den ! conductivity of snow 426 REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice 427 REAL, PARAMETER :: sea_cap=3995. ! specific heat capacity, water 428 REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice 429 430 ! control of snow and ice cover & freeze / melt (heights converted to kg/m2) 431 REAL :: snow_min!=sisno_min*sno_den !critical snow height 5 cm 432 REAL :: snow_wfact!=sisno_wfact ! max fraction of falling snow blown into ocean 433 REAL, PARAMETER :: ice_frac_min=0.005 434 REAL :: h_ice_min!=sithick_min ! min ice thickness 435 ! below ice_thin, priority is melt lateral / grow height 436 ! ice_thin is also height of new ice 437 REAL, PARAMETER :: h_ice_max=7 ! max ice height 438 ! Ice thickness parameter for lateral growth 439 REAL, PARAMETER :: h_ice_thick=1.5 440 REAL, PARAMETER :: h_ice_thin=0.15 441 442 ! albedo and radiation parameters 443 INTEGER, SAVE :: iflag_sic_albedo 444 ! albedo old or NEMO 445 REAL :: alb_sno_dry!=rn_alb_sdry !dry snow albedo 446 REAL :: alb_sno_wet!=rn_alb_smlt !wet snow albedo 447 REAL :: alb_ice_dry!=rn_alb_idry !dry thick ice 448 REAL :: alb_ice_wet!=rn_alb_imlt !melting thick ice 449 ! new (Toyoda 2020) albedo 450 ! Values for snow / ice, dry / melting, visible / near IR 451 REAL, PARAMETER :: alb_sdry_vis=0.98 452 REAL, PARAMETER :: alb_smlt_vis=0.88 453 REAL, PARAMETER :: alb_sdry_nir=0.7 454 REAL, PARAMETER :: alb_smlt_nir=0.55 455 REAL, PARAMETER :: alb_idry_vis=0.78 456 REAL, PARAMETER :: alb_imlt_vis=0.705 457 REAL, PARAMETER :: alb_idry_nir=0.36 458 REAL, PARAMETER :: alb_imlt_nir=0.285 459 REAL, PARAMETER :: h_ice_alb=0.5*ice_den ! height for full ice albedo 460 REAL, PARAMETER :: h_sno_alb=0.02*300 ! height for control of snow fraction 461 462 REAL :: pen_frac !=si_pen_frac !fraction of shortwave penetrating into the 463 ! ice (not snow). Should be visible only, not NIR 464 REAL :: pen_ext !=si_pen_ext !extinction length of penetrating shortwave (m-1) 465 466 ! HF from ocean below ice 467 ! REAL, PARAMETER :: fseaN=2.0 ! NH 468 ! REAL, PARAMETER :: fseaS=4.0 ! SH 469 !GG 344 470 345 471 #ifdef ISO … … 371 497 tsurf_tmp(:) = tsurf_in(:) 372 498 499 !GG 500 IF (iflag_seaice==0) THEN ! Old LMDZ sea ice surface 501 !GG 373 502 ! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal 374 503 CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd) … … 387 516 WHERE (snow > 0.0) cal = RCPD * calsno 388 517 ENDIF 518 519 !GG 520 ELSEIF (iflag_seaice==2) THEN 521 522 sno_den=sisno_den !mean snow density, kg/m3 523 ice_cond=sice_cond*ice_den !conductivity of ice 524 sno_cond=sisno_cond*sno_den ! conductivity of snow 525 snow_min=sisno_min*sno_den !critical snow height 5 cm 526 snow_wfact=sisno_wfact ! max fraction of falling snow blown into ocean 527 h_ice_min=sithick_min ! min ice thickness 528 alb_sno_dry=rn_alb_sdry !dry snow albedo 529 alb_sno_wet=rn_alb_smlt !wet snow albedo 530 alb_ice_dry=rn_alb_idry !dry thick ice 531 alb_ice_wet=rn_alb_imlt !melting thick ice 532 pen_frac=si_pen_frac !fraction of shortwave penetrating into the 533 pen_ext=si_pen_ext !extinction length of penetrating shortwave (m-1) 534 535 bilg(:)=0. 536 dif_grnd(:)=0. 537 beta(:) = 1. 538 fsic(:) = pctsrf(:,is_sic) 539 cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour 540 541 ! Surface, snow-ice and ice-ocean fluxes. 542 ! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd) 543 544 !write(*,*) 'radsol 1',radsol(1:100) 545 DO i=1,knon 546 ki=knindex(i) 547 IF (snow(i).GT.snow_min) THEN 548 ! 1 / snow-layer heat capacity 549 cal(i)=2.*RCPD/(snow(i)*ice_cap) 550 ! adjustment time-scale of conductive flux 551 dif_grnd(i) = cal(i) * sno_cond / snow(i) / RCPD 552 ! for conductive flux 553 f_cond_s = sno_cond * (tice(ki)-t_freeze) / snow(i) 554 radsol(i) = radsol(i)+f_cond_s 555 ! all shortwave flux absorbed 556 f_swpen=0. 557 ELSE ! bare ice. 558 f_cond_s = 0. 559 ! 1 / ice-layer heat capacity 560 cal(i) = 2.*RCPD/(hice(ki)*ice_den*ice_cap) 561 ! adjustment time-scale of conductive flux 562 dif_grnd(i) = cal(i) * ice_cond / (hice(ki)*ice_den) / RCPD 563 ! penetrative shortwave flux... 564 f_swpen=swnet(i)*pen_frac*exp(-pen_ext*hice(ki)) 565 radsol(i) = radsol(i)-f_swpen 566 ! GG no conductive flux in this case? 567 END IF 568 bilg(ki)=f_swpen-f_cond_s 569 END DO 570 571 endif 389 572 390 573 ! beta = 1.0 … … 423 606 ! 424 607 !**************************************************************************************** 608 !GG 609 if (iflag_seaice==0) then 610 !GG 425 611 #ifdef ISO 426 612 ! verif … … 502 688 alb2_new(:) = alb1_new(:) 503 689 690 !GG 691 else 692 693 DO i=1,knon 694 ki=knindex(i) 695 696 ! ocean-ice heat flux 697 fsea=fseaS 698 amax=amax_s 699 if (latitude(ki)>0) THEN 700 fsea=fseaN 701 amax=amax_n 702 ENDIF 703 704 IF (snow(i).GT.snow_min) THEN 705 ! snow conductive flux after calcul_fluxs (pos up) 706 f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i) 707 ! 1 / heat capacity and conductive timescale 708 uscap = 2. / ice_cap / (snow(i)+hice(ki)*ice_den) 709 ustau = uscap * ice_cond / (hice(ki)*ice_den) 710 ! update ice temp 711 tice(ki) = (tice(ki) + dtime*(ustau*t_freeze - uscap*f_cond_s)) & 712 / (1 + dtime*ustau) 713 ELSE ! bare ice 714 tice(ki)=tsurf_new(i) 715 ENDIF 716 ! ice conductive flux (pos up) 717 f_cond_i = ice_cond * (t_freeze-tice(ki)) / (hice(ki)*ice_den) 718 f_bot(i) = fsea - f_cond_i 719 fcdi(ki) = f_cond_i - fsea 720 fcds(i) = f_cond_s 721 !bilg(ki) = bilg(ki)+f_cond_i 722 END DO 723 724 !**************************************************************************************** 725 ! 2) Update snow and ice surface : thickness 726 !**************************************************************************************** 727 728 IF (iflag_seaice==1) THEN 729 ! Read from limit 730 CALL limit_read_hice(knon,knindex,hice) 731 ENDIF 732 ! Formula Krinner et al. 1997 : h = (0.2 + 3.8(f_min**2))(1 + 2(f- f_min)) 733 734 735 736 DO i=1,knon 737 ki=knindex(i) 738 IF (precip_snow(i) > 0.) THEN 739 snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki))) 740 END IF 741 ! snow and ice sublimation 742 IF (evap(i) > 0.) THEN 743 snow_evap = MIN (snow(i) / dtime, evap(i)) 744 snow(i) = snow(i) - snow_evap * dtime 745 snow(i) = MAX(0.0, snow(i)) 746 IF (iflag_seaice==2) THEN 747 hice(ki) = MAX(0.0,hice(ki)-(evap(i)-snow_evap)*dtime/ice_den) 748 ENDIF 749 ENDIF 750 ! Melt / Freeze snow from above if Tsurf>0 751 IF (tsurf_new(i).GT.t_melt) THEN 752 ! energy available for melting snow (in kg of melted snow /m2) 753 e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. & 754 /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i)) 755 ! remove snow 756 tice_i=tice(ki) 757 IF (snow(i).GT.e_melt) THEN 758 snow(i)=snow(i)-e_melt 759 tsurf_new(i)=t_melt 760 ELSE ! all snow is melted 761 ! add remaining heat flux to ice 762 e_melt=e_melt-snow(i) 763 tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*hice(ki)*ice_den) 764 tsurf_new(i)=tice(ki) 765 END IF 766 dtice_melt(ki)=(tice(ki)-tice_i)/dtime 767 END IF 768 ! Bottom melt / grow 769 ! bottom freeze if bottom flux (cond + oce-ice) <0 770 IF (iflag_seaice==2) THEN 771 IF (f_bot(i).LT.0) THEN 772 ! larger fraction of bottom growth 773 frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thick) & 774 / (h_ice_max-h_ice_thick))) 775 ! quantity of new ice (formed at mean ice temp) 776 e_melt= -f_bot(i) * dtime * fsic(ki) & 777 / (ice_lat+ice_cap/2.*(t_freeze-tice(ki))) 778 ! first increase height to h_thick 779 dhsic=MAX(0.,MIN(h_ice_thick-hice(ki),e_melt/(fsic(ki)*ice_den))) 780 hice_i=hice(ki) 781 hice(ki)=dhsic+hice(ki) 782 e_melt=e_melt-fsic(ki)*dhsic 783 IF (e_melt.GT.0.) THEN 784 ! frac_mf fraction used for lateral increase 785 dfsic=MIN(amax-fsic(ki),e_melt*frac_mf/ (hice(ki)*ice_den) ) 786 ! No lateral growth -> forced ocean 787 !fsic(ki)=fsic(ki)+dfsic 788 e_melt=e_melt-dfsic*hice(ki)*ice_den 789 ! rest used to increase height 790 hice(ki)=MIN(h_ice_max,hice(ki)+e_melt/( fsic(ki) * ice_den ) ) 791 END IF 792 dh_basal_growth(ki)=(hice(ki)-hice_i)/dtime 793 794 ! melt from below if bottom flux >0 795 ELSE 796 ! larger fraction of lateral melt from warm ocean 797 frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thin) & 798 / (h_ice_thick-h_ice_thin))) 799 ! bring ice to freezing and melt from below 800 ! quantity of melted ice 801 e_melt= f_bot(i) * dtime * fsic(ki) & 802 / (ice_lat+ice_cap/2.*(tice(ki)-t_freeze)) 803 ! first decrease height to h_thick 804 hice_i=hice(ki) 805 dhsic=MAX(0.,MIN(hice(ki)-h_ice_thick,e_melt/(fsic(ki)*ice_den))) 806 hice(ki)=hice(ki)-dhsic 807 e_melt=e_melt-fsic(ki)*dhsic*ice_den 808 809 IF (e_melt.GT.0) THEN 810 ! frac_mf fraction used for height decrease 811 dhsic=MAX(0.,MIN(hice(ki)-h_ice_min,e_melt/ice_den*frac_mf/fsic(ki))) 812 hice(ki)=hice(ki)-dhsic 813 e_melt=e_melt-fsic(ki)*dhsic*ice_den 814 ! rest used to decrease fraction (up to 0!) 815 dfsic=MIN(fsic(ki),e_melt/(hice(ki)*ice_den)) 816 ! Remaining heat not used if everything melted 817 e_melt=e_melt-dfsic*hice(ki)*ice_den 818 ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime 819 END IF 820 dh_basal_melt(ki)=(hice(ki)-hice_i)/dtime 821 END IF 822 END IF 823 824 ! melt ice from above if Tice>0 825 tice_i=tice(ki) 826 IF (tice(ki).GT.t_melt) THEN 827 IF (iflag_seaice==2) THEN 828 ! quantity of ice melted (kg/m2) 829 e_melt=MAX(hice(ki)*ice_den*(tice(ki)-t_melt)*ice_cap/2. & 830 /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0) 831 ! melt from above, height only 832 hice_i=hice(ki) 833 dhsic=MIN(hice(ki)-h_ice_min,e_melt/ice_den) 834 dh_top_melt(i)=dhsic 835 e_melt=e_melt-dhsic 836 IF (e_melt.GT.0) THEN 837 ! lateral melt if ice too thin 838 dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/(h_ice_min*ice_den)*fsic(ki)) 839 ! if all melted do nothing with remaining heat 840 e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min*ice_den) 841 ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime 842 END IF 843 hice(ki)=hice(ki)-dhsic 844 dh_top_melt(ki)=(hice(ki)-hice_i)/dtime 845 ! surface temperature at melting point 846 END IF 847 tice(ki)=t_melt 848 tsurf_new(i)=t_melt 849 END IF 850 dtice_melt(ki)=dtice_melt(ki)+tice(ki)-tice_i 851 852 ! convert snow to ice if below floating line 853 h_test=(hice(ki)*ice_den+snow(i))-hice(ki)*sea_den 854 IF ((h_test.GT.0.).AND.(hice(ki).GT.h_ice_min)) THEN !snow under water 855 ! extra snow converted to ice (with added frozen sea water) 856 IF (iflag_seaice==2) THEN 857 dhsic=h_test/(sea_den-ice_den+sno_den) 858 hice(ki)=hice(ki)+dhsic 859 ENDIF 860 snow(i)=snow(i)-dhsic*sno_den 861 ! available energy (freeze sea water + bring to tice) 862 e_melt=dhsic*ice_den*(1.-sno_den/ice_den)*(ice_lat+ & 863 ice_cap/2.*(t_freeze-tice(ki))) 864 ! update ice temperature 865 tice_i=tice(ki) 866 tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+hice(ki)*ice_den) 867 IF (iflag_seaice==2) THEN 868 dh_snow2sic(ki)=dhsic/dtime 869 END IF 870 dtice_snow2sic(ki)=(tice(ki)-tice_i)/dtime 871 END IF 872 END DO 873 874 !write(*,*) 'hice 2',hice(1:100) 875 !write(*,*) 'tice 2',tice(1:100) 876 877 iflag_sic_albedo=iflag_seaice_alb 878 879 !******************************************************************************* 880 ! 3) cumulate ice-ocean fluxes, update tslab, lateral grow 881 !***********************************************o******************************* 882 !cumul fluxes. 883 bilg_cumul(:)=bilg_cumul(:)+bilg(:)/float(cpl_pas) 884 IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab 885 bilg_cumul(:)=0. 886 END IF ! coupling time 887 888 ! write(*,*) 'hice 3',hice(1:100) 889 ! write(*,*) 'tice 3',tice(1:100) 890 !tests ice fraction 891 WHERE (fsic.LT.ice_frac_min) 892 tice=t_melt 893 hice=h_ice_min 894 END WHERE 895 896 !write(*,*) 'hice 4',hice(1:100) 897 !write(*,*) 'tice 4',tice(1:100) 898 899 endif 900 901 !**************************************************************************************** 902 ! 4) Compute sea-ice and snow albedo 903 !**************************************************************************************** 904 IF (iflag_seaice > 0) THEN 905 SELECT CASE (iflag_sic_albedo) 906 CASE(0) 907 ! old slab albedo : single value. age of snow, melt ponds. 908 DO i=1,knon 909 ki=knindex(i) 910 ! snow albedo: update snow age 911 IF (snow(i).GT.0.0001) THEN 912 agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)& 913 * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.) 914 ELSE 915 agesno(i)=0.0 916 END IF 917 ! snow albedo 918 alb_snow=alb_sno_wet+(alb_sno_dry-alb_sno_wet)*EXP(-agesno(i)/50.) 919 ! ice albedo (varies with ice tkickness and temp) 920 alb_ice=MAX(0.0,0.13*LOG(100.*hice(ki))+0.1) 921 !alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1) 922 IF (tice(ki).GT.t_freeze-0.01) THEN 923 alb_ice=MIN(alb_ice,alb_ice_wet) 924 ELSE 925 alb_ice=MIN(alb_ice,alb_ice_dry) 926 END IF 927 ! pond albedo 928 alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) 929 ! pond fraction 930 frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) 931 ! snow fraction 932 frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min)) 933 ! ice fraction 934 frac_ice=MAX(0.0,1.-frac_pond-frac_snow) 935 ! total albedo 936 alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond 937 END DO 938 alb2_new(:) = alb1_new(:) 939 940 CASE(1) 941 ! New visible and IR albedos, dry / melting snow 942 ! based on Toyoda et al, 2020 943 DO i=1,knon 944 ki=knindex(i) 945 ! snow fraction 946 frac_snow = snow(i) / (snow(i) + h_sno_alb) 947 ! dependence of ice albedo with ice thickness 948 frac_ice = MIN(1.,ATAN(4.*hice(ki)*ice_den) / ATAN(4.*h_ice_alb)) 949 ! Total (for ice, min = 0.066 = alb_ocean) 950 IF (tice(ki).GT.t_melt) THEN 951 alb_ice = 0.066 + (alb_imlt_vis - 0.066)*frac_ice 952 alb1_new(i)=alb_smlt_vis*frac_snow + alb_ice*(1.-frac_snow) 953 alb_ice = 0.066 + (alb_imlt_nir - 0.066)*frac_ice 954 alb2_new(i)=alb_smlt_nir*frac_snow + alb_ice*(1.-frac_snow) 955 ELSEIF (tice(ki).GT.t_melt - 1.) THEN 956 frac_pond = tice(ki) - t_freeze 957 alb_snow = alb_smlt_vis*frac_pond + alb_sdry_vis*(1.-frac_pond) 958 alb_ice = alb_imlt_vis*frac_pond + alb_idry_vis*(1.-frac_pond) 959 alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice 960 alb1_new(i)= alb_snow*frac_snow + alb_ice*(1.-frac_snow) 961 alb_snow = alb_smlt_nir*frac_pond + alb_sdry_nir*(1.-frac_pond) 962 alb_ice = alb_imlt_nir*frac_pond + alb_idry_nir*(1.-frac_pond) 963 alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice 964 alb2_new(i)= alb_snow*frac_snow + alb_ice*(1.-frac_snow) 965 ELSE 966 alb_ice = 0.066 + (alb_idry_vis - 0.066)*frac_ice 967 alb1_new(i)=alb_sdry_vis*frac_snow + alb_ice*(1.-frac_snow) 968 alb_ice = 0.066 + (alb_idry_nir - 0.066)*frac_ice 969 alb2_new(i)=alb_sdry_nir*frac_snow + alb_ice*(1.-frac_snow) 970 ENDIF 971 END DO 972 973 CASE(2) 974 ! LIM3 scheme. Uses clear sky / overcast value, with 50% clear sky 975 z1_i = 1.5 * ice_den 976 z2_i = 0.05 * ice_den 977 zlog = 1. / (LOG(1.5) - LOG(0.05)) 978 z1_s = 1. / (0.025 * sno_den) 979 DO i=1,knon 980 ki=knindex(i) 981 ! temperature above / below 0 982 IF (tice(ki).GE.t_melt) THEN 983 alb_ice = alb_ice_wet 984 alb_snow = alb_sno_wet 985 ELSE 986 alb_ice = alb_ice_dry 987 alb_snow = alb_sno_dry 988 ENDIF 989 ! ice thickness 990 IF (hice(ki)*ice_den.LT.z2_i) THEN 991 alb_ice = 0.066 + 0.114 * hice(ki)*ice_den / z2_i 992 ELSEIF (hice(ki)*ice_den.LT.z1_i) THEN 993 alb_ice = alb_ice + (0.18 - alb_ice) & 994 * (LOG(z1_i) - LOG(hice(ki)*ice_den)) * zlog 995 ENDIF 996 ! ice or snow depending on snow thickness 997 alb_snow = alb_snow - (alb_snow -alb_ice) * EXP(- snow(i) * z1_s) 998 ! Effect of clouds (polynomial fit with 50% clouds) 999 alb1_new(i) = alb_snow - 0.5 * (-0.1010 * alb_snow*alb_snow & 1000 + 0.1933*alb_snow - 0.0148) 1001 alb2_new(i) = alb1_new(i) 1002 END DO 1003 1004 CASE(3) 1005 CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) 1006 WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0. 1007 alb1_new(:) = 0.0 1008 DO i=1, knon 1009 zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0))) 1010 alb1_new(i) = alb_neig(i) * zfra + 0.6 * (1.0-zfra) 1011 ENDDO 1012 alb2_new(:) = alb1_new(:) 1013 1014 print*,'alb_neig=',alb_neig 1015 print*,'zfra=',zfra 1016 print*,'snow=',snow 1017 print*,'alb1_new=',alb1_new 1018 print*,'alb2_new=',alb2_new 1019 END SELECT 1020 END IF 1021 ! ------ End Albedo ---------- 1022 1023 !GG 504 1024 END SUBROUTINE ocean_forced_ice 505 1025
Note: See TracChangeset
for help on using the changeset viewer.
