Changeset 12937 for NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE
- Timestamp:
- 2020-05-15T18:15:25+02:00 (4 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE
- Files:
-
- 22 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/ice.F90
r11715 r12937 70 70 !! a_ip | - | Ice pond concentration | | 71 71 !! v_ip | - | Ice pond volume per unit area| m | 72 !! v_il | v_il_1d | Ice pond lid volume per area | m | 72 73 !! | 73 74 !!-------------|-------------|---------------------------------|-------| … … 85 86 !! t_su ! t_su_1d | Sea ice surface temperature ! K | 86 87 !! h_ip | h_ip_1d | Ice pond thickness | m | 88 !! h_il | h_il_1d | Ice pond lid thickness | m | 87 89 !! | 88 90 !! notes: the ice model only sees a bulk (i.e., vertically averaged) | … … 112 114 !! hm_ip | - | Mean ice pond depth | m | 113 115 !! vt_ip | - | Total ice pond vol. per unit area| m | 116 !! hm_il | - | Mean ice pond lid depth | m | 117 !! vt_il | - | Total ice pond lid vol. per area | m | 114 118 !!===================================================================== 115 119 … … 151 155 INTEGER , PUBLIC :: nn_nevp !: number of iterations for subcycling 152 156 REAL(wp), PUBLIC :: rn_relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 157 LOGICAL , PUBLIC :: ln_rhg_chkcvg !: check ice rheology convergence 153 158 ! 154 159 ! !!** ice-advection namelist (namdyn_adv) ** … … 190 195 ! !!** ice-ponds namelist (namthd_pnd) 191 196 LOGICAL , PUBLIC :: ln_pnd !: Melt ponds (T) or not (F) 192 LOGICAL , PUBLIC :: ln_pnd_H12 !: Melt ponds scheme from Holland et al 2012 197 LOGICAL , PUBLIC :: ln_pnd_H12 !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010) 198 LOGICAL, PUBLIC :: ln_pnd_lids !: Allow ponds to have frozen lids 199 LOGICAL, PUBLIC :: ln_pnd_flush !: Allow ponds to flush thru the ice 200 REAL(wp), PUBLIC :: rn_apnd_min !: Minimum ice fraction that contributes to melt ponds 201 REAL(wp), PUBLIC :: rn_apnd_max !: Maximum ice fraction that contributes to melt ponds 193 202 LOGICAL , PUBLIC :: ln_pnd_CST !: Melt ponds scheme with constant fraction and depth 194 203 REAL(wp), PUBLIC :: rn_apnd !: prescribed pond fraction (0<rn_apnd<1) … … 271 280 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_snw !: heat flux for snow melt [W.m-2] 272 281 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2] 273 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_err_rem !: heat flux error after heat remapping => must be 0 [W.m-2]274 282 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qt_atm_oi !: heat flux at the interface atm-[oce+ice] [W.m-2] 275 283 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qt_oce_ai !: heat flux at the interface oce-[atm+ice] [W.m-2] … … 331 339 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_ip !: melt pond volume per grid cell area [m] 332 340 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip_frac !: melt pond fraction (a_ip/a_i) 341 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip_eff !: melt pond effective fraction (not covered up by lid) (a_ip/a_i) 333 342 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_ip !: melt pond depth [m] 343 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_il !: melt pond lid volume [m] 344 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_il !: melt pond lid thickness [m] 334 345 335 346 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_ip !: total melt pond concentration 336 347 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hm_ip !: mean melt pond depth [m] 337 348 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_ip !: total melt pond volume per gridcell area [m] 349 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hm_il !: mean melt pond lid depth [m] 350 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_il !: total melt pond lid volume per gridcell area [m] 338 351 339 352 !!---------------------------------------------------------------------- 340 353 !! * Old values of global variables 341 354 !!---------------------------------------------------------------------- 342 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s_b, v_i_b, h_s_b, h_i_b , h_ip_b!: snow and ice volumes/thickness343 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_b, sv_i_b , oa_i_b!:344 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s_b 345 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i_b 346 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice_b, v_ice_b 347 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i_b 355 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness 356 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_b, sv_i_b !: 357 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s_b !: snow heat content 358 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i_b !: ice temperatures 359 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice_b, v_ice_b !: ice velocity 360 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i_b !: ice concentration (total) 348 361 349 362 !!---------------------------------------------------------------------- … … 386 399 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice_top !: Surface conduction flux (W/m2) 387 400 401 !!---------------------------------------------------------------------- 402 !! * Only for atmospheric coupling 403 !!---------------------------------------------------------------------- 404 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Ice fractional area at last coupling time 388 405 ! 389 406 !!---------------------------------------------------------------------- … … 400 417 INTEGER :: ice_alloc 401 418 ! 402 INTEGER :: ierr(1 6), ii419 INTEGER :: ierr(17), ii 403 420 !!----------------------------------------------------------------- 404 421 ierr(:) = 0 … … 424 441 & hfx_sum (jpi,jpj) , hfx_bom (jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , & 425 442 & hfx_opw (jpi,jpj) , hfx_thd (jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) , & 426 & hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj) ,STAT=ierr(ii) )443 & hfx_err_dif(jpi,jpj) , wfx_err_sub(jpi,jpj) , STAT=ierr(ii) ) 427 444 428 445 ! * Ice global state variables … … 448 465 449 466 ii = ii + 1 450 ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) ) 451 452 ii = ii + 1 453 ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) ) 467 ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), & 468 & v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , STAT = ierr(ii) ) 469 470 ii = ii + 1 471 ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , hm_il(jpi,jpj) , vt_il(jpi,jpj) , STAT = ierr(ii) ) 454 472 455 473 ! * Old values of global variables 456 474 ii = ii + 1 457 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl) , h_i_b(jpi,jpj,jpl), h_ip_b(jpi,jpj,jpl),&458 & a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , 459 & oa_i_b(jpi,jpj,jpl) ,STAT=ierr(ii) )475 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl) , h_i_b(jpi,jpj,jpl), & 476 & a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 477 & STAT=ierr(ii) ) 460 478 461 479 ii = ii + 1 … … 481 499 ALLOCATE( t_si(jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj,jpl) , qcn_ice_top(jpi,jpj,jpl) , STAT = ierr(ii) ) 482 500 501 ! * For atmospheric coupling 502 ii = ii + 1 503 ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(ii) ) 504 483 505 ice_alloc = MAXVAL( ierr(:) ) 484 506 IF( ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' ) 485 507 ! 508 486 509 END FUNCTION ice_alloc 487 510 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/ice1d.F90
r11715 r12937 51 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_snw_1d 52 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_dyn_1d 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_err_rem_1d54 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_err_dif_1d 55 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qt_oce_ai_1d … … 124 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: oa_i_1d !: 125 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: o_i_1d !: 126 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_1d !: 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_1d !: ice ponds 127 126 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_ip_1d !: 128 127 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_ip_1d !: 129 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_frac_1d !: 128 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_il_1d !: Ice pond lid 129 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_il_1d !: 130 130 131 131 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s … … 157 157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_ip_2d 158 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_ip_2d 159 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_il_2d 159 160 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_su_2d 160 161 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_2d … … 189 190 & hfx_thd_1d(jpij) , hfx_spr_1d (jpij) , & 190 191 & hfx_snw_1d(jpij) , hfx_sub_1d (jpij) , & 191 & hfx_res_1d(jpij) , hfx_err_ rem_1d(jpij) , hfx_err_dif_1d(jpij) , qt_oce_ai_1d(jpij), STAT=ierr(ii) )192 & hfx_res_1d(jpij) , hfx_err_dif_1d(jpij) , qt_oce_ai_1d(jpij), STAT=ierr(ii) ) 192 193 ! 193 194 ii = ii + 1 … … 208 209 & dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) , & 209 210 & dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d (jpij) , s_i_new (jpij) , & 210 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , 211 & h_i p_1d (jpij) , a_ip_frac_1d(jpij) ,&211 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , v_il_1d (jpij) , & 212 & h_il_1d (jpij) , h_ip_1d (jpij) , & 212 213 & sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d (jpij) , STAT=ierr(ii) ) 213 214 ! … … 226 227 ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) , & 227 228 & v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) , & 228 & a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , 229 & a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , v_il_2d(jpij,jpl) , & 229 230 & STAT=ierr(ii) ) 230 231 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icealb.F90
r11715 r12937 45 45 CONTAINS 46 46 47 SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, p alb_cs, palb_os)47 SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, pcloud_fra, palb_ice ) 48 48 !!---------------------------------------------------------------------- 49 49 !! *** ROUTINE ice_alb *** … … 97 97 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pafrac_pnd ! melt pond relative fraction (per unit ice area) 98 98 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: ph_pnd ! melt pond depth 99 REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_cs ! albedo of ice under clear sky100 REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_ os ! albedo of ice under overcast sky99 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pcloud_fra ! cloud fraction 100 REAL(wp), INTENT( out), DIMENSION(:,:,:) :: palb_ice ! albedo of ice 101 101 ! 102 102 INTEGER :: ji, jj, jl ! dummy loop indices … … 106 106 REAL(wp) :: zalb_ice, zafrac_ice ! bare sea ice albedo & relative ice fraction 107 107 REAL(wp) :: zalb_snw, zafrac_snw ! snow-covered sea ice albedo & relative snow fraction 108 REAL(wp) :: zalb_cs, zalb_os ! albedo of ice under clear/overcast sky 108 109 !!--------------------------------------------------------------------- 109 110 ! … … 119 120 DO jj = 1, jpj 120 121 DO ji = 1, jpi 121 ! !--- Specific snow, ice and pond fractions (for now, we prevent melt ponds and snow at the same time) 122 IF( ph_snw(ji,jj,jl) == 0._wp ) THEN 122 !---------------------------------------------! 123 !--- Specific snow, ice and pond fractions ---! 124 !---------------------------------------------! 125 IF( ph_snw(ji,jj,jl) == 0._wp ) THEN !--- no snow : we prevent melt ponds and snow at the same time (for now) 123 126 zafrac_snw = 0._wp 124 127 IF( ld_pnd_alb ) THEN … … 129 132 zafrac_ice = 1._wp - zafrac_pnd 130 133 ELSE 131 zafrac_snw = 1._wp ! Snowfully "shades" melt ponds and ice134 zafrac_snw = 1._wp !--- snow : fully "shades" melt ponds and ice 132 135 zafrac_pnd = 0._wp 133 136 zafrac_ice = 0._wp 134 137 ENDIF 135 138 ! 139 !---------------! 140 !--- Albedos ---! 141 !---------------! 136 142 ! !--- Bare ice albedo (for hi > 150cm) 137 143 IF( ld_pnd_alb ) THEN … … 161 167 ENDIF 162 168 ! !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions 163 palb_os(ji,jj,jl) = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1) 164 ! 165 palb_cs(ji,jj,jl) = palb_os(ji,jj,jl) & 166 & - ( - 0.1010 * palb_os(ji,jj,jl) * palb_os(ji,jj,jl) & 167 & + 0.1933 * palb_os(ji,jj,jl) - 0.0148 ) * tmask(ji,jj,1) 168 ! 169 zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1) 170 ! 171 zalb_cs = zalb_os - ( - 0.1010 * zalb_os * zalb_os & 172 & + 0.1933 * zalb_os - 0.0148 ) * tmask(ji,jj,1) 173 ! 174 ! albedo depends on cloud fraction because of non-linear spectral effects 175 palb_ice(ji,jj,jl) = ( 1._wp - pcloud_fra(ji,jj) ) * zalb_cs + pcloud_fra(ji,jj) * zalb_os 176 169 177 END DO 170 178 END DO -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icectl.F90
r11715 r12937 625 625 WRITE(numout,*) ' e_snow : ', e_s(ji,jj,1,jl) , ' e_snow_b : ', e_s_b(ji,jj,1,jl) 626 626 WRITE(numout,*) ' sv_i : ', sv_i(ji,jj,jl) , ' sv_i_b : ', sv_i_b(ji,jj,jl) 627 WRITE(numout,*) ' oa_i : ', oa_i(ji,jj,jl) , ' oa_i_b : ', oa_i_b(ji,jj,jl)628 627 END DO !jl 629 628 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn.F90
r11715 r12937 99 99 WHERE( a_ip(:,:,:) >= epsi20 ) 100 100 h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 101 h_il(:,:,:) = v_il(:,:,:) / a_ip(:,:,:) 101 102 ELSEWHERE 102 103 h_ip(:,:,:) = 0._wp 104 h_il(:,:,:) = 0._wp 103 105 END WHERE 104 106 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn_adv.F90
r11715 r12937 84 84 ! !-----------------------! 85 85 CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 86 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )86 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, v_il, e_s, e_i ) 87 87 ! !-----------------------! 88 88 CASE( np_advPRA ) ! PRATHER scheme ! 89 89 ! !-----------------------! 90 90 CALL ice_dyn_adv_pra( kt, u_ice, v_ice, & 91 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )91 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, v_il, e_s, e_i ) 92 92 END SELECT 93 93 -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn_adv_pra.F90
r11715 r12937 44 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxap , syap , sxxap , syyap , sxyap ! melt pond fraction 45 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvp , syvp , sxxvp , syyvp , sxyvp ! melt pond volume 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvl , syvl , sxxvl , syyvl , sxyvl ! melt pond lid volume 46 47 47 48 !! * Substitutions … … 55 56 56 57 SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, & 57 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )58 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 58 59 !!---------------------------------------------------------------------- 59 60 !! ** routine ice_dyn_adv_pra ** … … 78 79 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 79 80 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 81 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid thickness 80 82 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 81 83 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 88 90 REAL(wp), DIMENSION(jpi,jpj,1) :: z0opw 89 91 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi 90 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp 92 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp, z0vl 91 93 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 92 94 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei … … 129 131 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 130 132 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 133 IF ( ln_pnd_lids ) THEN 134 z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:) ! Melt pond lid volume 135 ENDIF 131 136 ENDIF 132 137 END DO … … 167 172 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 168 173 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 174 IF ( ln_pnd_lids ) THEN 175 CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 176 CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 177 ENDIF 169 178 ENDIF 170 179 END DO … … 202 211 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 203 212 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 213 IF ( ln_pnd_lids ) THEN 214 CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 215 CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 216 ENDIF 204 217 ENDIF 205 218 END DO … … 225 238 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 226 239 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 240 IF ( ln_pnd_lids ) THEN 241 pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 242 ENDIF 227 243 ENDIF 228 244 END DO … … 231 247 ! Remove negative values (conservation is ensured) 232 248 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 233 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )249 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 234 250 ! 235 251 ! --- Ensure snow load is not too big --- ! … … 651 667 & sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) , & 652 668 & sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) , & 653 & sxap(jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , & 654 & sxvp(jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , & 669 & sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , & 670 & sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , & 671 & sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) , & 655 672 ! 656 673 & sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & … … 765 782 CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 766 783 CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 784 ! 785 IF ( ln_pnd_lids ) THEN ! melt pond lid volume 786 CALL iom_get( numrir, jpdom_autoglo, 'sxvl' , sxvl ) 787 CALL iom_get( numrir, jpdom_autoglo, 'syvl' , syvl ) 788 CALL iom_get( numrir, jpdom_autoglo, 'sxxvl', sxxvl ) 789 CALL iom_get( numrir, jpdom_autoglo, 'syyvl', syyvl ) 790 CALL iom_get( numrir, jpdom_autoglo, 'sxyvl', sxyvl ) 791 ENDIF 767 792 ENDIF 768 793 ! … … 780 805 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content 781 806 IF( ln_pnd_H12 ) THEN 782 sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction 783 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume 807 sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction 808 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume 809 IF ( ln_pnd_lids ) THEN 810 sxvl = 0._wp; syvl = 0._wp ; sxxvl = 0._wp ; syyvl = 0._wp ; sxyvl = 0._wp ! melt pond lid volume 811 ENDIF 784 812 ENDIF 785 813 ENDIF … … 862 890 CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 863 891 CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp ) 892 ! 893 IF ( ln_pnd_lids ) THEN ! melt pond lid volume 894 CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl ) 895 CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl ) 896 CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl ) 897 CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl ) 898 CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl ) 899 ENDIF 864 900 ENDIF 865 901 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn_adv_umx.F90
r11715 r12937 60 60 61 61 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 62 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )62 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 63 63 !!---------------------------------------------------------------------- 64 64 !! *** ROUTINE ice_dyn_adv_umx *** … … 85 85 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond concentration 86 86 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 87 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid volume 87 88 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 88 89 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 334 335 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 335 336 & zhvar, pv_ip, zua_ups, zva_ups ) 337 ! lid 338 IF ( ln_pnd_lids ) THEN 339 zamsk = 0._wp 340 zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:) 341 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 342 & zhvar, pv_il, zua_ups, zva_ups ) 343 ENDIF 336 344 ENDIF 337 345 ! … … 350 358 ! Remove negative values (conservation is ensured) 351 359 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 352 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )360 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 353 361 ! 354 362 ! Make sure ice thickness is not too big -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn_rdgrft.F90
r11715 r12937 503 503 REAL(wp) :: airdg1, oirdg1, aprdg1, virdg1, sirdg1 504 504 REAL(wp) :: airft1, oirft1, aprft1 505 REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg ! area etc of new ridges506 REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft ! area etc of rafted ice505 REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg ! area etc of new ridges 506 REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft ! area etc of rafted ice 507 507 ! 508 508 REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges … … 581 581 aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 582 582 vprft (ji) = v_ip_2d(ji,jl1) * afrft 583 IF ( ln_pnd_lids ) THEN 584 vlrdg (ji) = v_il_2d(ji,jl1) * afrdg 585 vlrft (ji) = v_il_2d(ji,jl1) * afrft 586 ENDIF 583 587 ENDIF 584 588 … … 610 614 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1 611 615 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 616 IF ( ln_pnd_lids ) THEN 617 v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - vlrdg(ji) - vlrft(ji) 618 ENDIF 612 619 ENDIF 613 620 ENDIF … … 706 713 a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea & 707 714 & + aprft2(ji) * rn_fpndrft * zswitch(ji) ) 715 IF ( ln_pnd_lids ) THEN 716 v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + ( vlrdg (ji) * rn_fpndrdg * fvol (ji) & 717 & + vlrft (ji) * rn_fpndrft * zswitch(ji) ) 718 ENDIF 708 719 ENDIF 709 720 … … 736 747 !---------------- 737 748 ! In case ridging/rafting lead to very small negative values (sometimes it happens) 738 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )749 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d ) 739 750 ! 740 751 END SUBROUTINE rdgrft_shift … … 848 859 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 849 860 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 861 CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) ) 850 862 DO jl = 1, jpl 851 863 DO jk = 1, nlay_s … … 874 886 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 875 887 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 888 CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) ) 876 889 DO jl = 1, jpl 877 890 DO jk = 1, nlay_s -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn_rhg.F90
r11715 r12937 110 110 INTEGER :: ios, ioptio ! Local integer output status for namelist read 111 111 !! 112 NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast 112 NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg 113 113 !!------------------------------------------------------------------- 114 114 ! … … 126 126 WRITE(numout,*) '~~~~~~~~~~~~~~~' 127 127 WRITE(numout,*) ' Namelist : namdyn_rhg:' 128 WRITE(numout,*) ' rheology EVP (icedyn_rhg_evp) ln_rhg_EVP = ', ln_rhg_EVP 129 WRITE(numout,*) ' use adaptive EVP (aEVP) ln_aEVP = ', ln_aEVP 130 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 131 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 132 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 133 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 128 WRITE(numout,*) ' rheology EVP (icedyn_rhg_evp) ln_rhg_EVP = ', ln_rhg_EVP 129 WRITE(numout,*) ' use adaptive EVP (aEVP) ln_aEVP = ', ln_aEVP 130 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 131 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 132 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 133 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 134 WRITE(numout,*) ' check convergence of rheology ln_rhg_chkcvg = ', ln_rhg_chkcvg 134 135 ENDIF 135 136 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn_rhg_evp.F90
r11715 r12937 41 41 USE prtctl ! Print control 42 42 43 USE netcdf ! NetCDF library for convergence test 43 44 IMPLICIT NONE 44 45 PRIVATE … … 49 50 !! * Substitutions 50 51 # include "vectopt_loop_substitute.h90" 52 53 !! for convergence tests 54 INTEGER :: ncvgid ! netcdf file id 55 INTEGER :: nvarid ! netcdf variable id 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 51 57 !!---------------------------------------------------------------------- 52 58 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 125 131 REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 126 132 ! 127 REAL(wp) :: zresm ! Maximal error on ice velocity128 133 REAL(wp) :: zintb, zintn ! dummy argument 129 134 REAL(wp) :: zfac_x, zfac_y … … 141 146 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 142 147 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 143 !!$ REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence144 148 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 145 149 ! ! ocean surface (ssh_m) if ice is not embedded … … 160 164 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 161 165 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 166 !! --- check convergence 167 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 162 168 !! --- diags 163 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00164 169 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 165 170 !! --- SIMIP diags … … 174 179 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 175 180 ! 176 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 181 ! for diagnostics and convergence tests 182 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 183 DO jj = 1, jpj 184 DO ji = 1, jpi 185 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 186 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 187 END DO 188 END DO 189 ! 190 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... 177 191 !------------------------------------------------------------------------------! 178 192 ! 0) mask at F points for the ice … … 345 359 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 346 360 ! 347 !!$ IF(ln_ctl) THEN ! Convergence test 348 !!$ DO jj = 1, jpjm1 349 !!$ zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 350 !!$ zv_ice(:,jj) = v_ice(:,jj) 351 !!$ END DO 352 !!$ ENDIF 361 ! convergence test 362 IF(ln_rhg_chkcvg) THEN 363 DO jj = 1, jpj 364 DO ji = 1, jpi 365 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 366 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 367 END DO 368 END DO 369 ENDIF 353 370 354 371 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! … … 667 684 ENDIF 668 685 669 !!$ IF(ln_ctl) THEN ! Convergence test 670 !!$ DO jj = 2 , jpjm1 671 !!$ zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 672 !!$ END DO 673 !!$ zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 674 !!$ CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 675 !!$ ENDIF 686 ! convergence test 687 IF(ln_rhg_chkcvg) THEN 688 CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 689 ENDIF 676 690 ! 677 691 ! ! ==================== ! … … 734 748 ! 5) diagnostics 735 749 !------------------------------------------------------------------------------! 736 DO jj = 1, jpj737 DO ji = 1, jpi738 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice739 END DO740 END DO741 742 750 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 743 751 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & … … 796 804 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 797 805 ENDIF 798 806 799 807 ! --- SIMIP --- ! 800 808 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & … … 852 860 ENDIF 853 861 ! 862 ! --- convergence tests --- ! 863 IF( ln_rhg_chkcvg ) THEN 864 IF( iom_use('uice_cvg') ) THEN ! output: u(t=nn_nevp) - u(t=nn_nevp-1) 865 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 866 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 867 ENDIF 868 ENDIF 869 ! 870 DEALLOCATE( zmsk00, zmsk15 ) 871 ! 854 872 END SUBROUTINE ice_dyn_rhg_evp 873 874 875 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 876 !!---------------------------------------------------------------------- 877 !! *** ROUTINE rhg_cvg *** 878 !! 879 !! ** Purpose : check convergence of oce rheology 880 !! 881 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity 882 !! during the sub timestepping of rheology so as: 883 !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 884 !! This routine is called every sub-iteration, so it is cpu expensive 885 !! 886 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 887 !!---------------------------------------------------------------------- 888 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 889 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 890 !! 891 INTEGER :: it, idtime, istatus 892 INTEGER :: ji, jj ! dummy loop indices 893 REAL(wp) :: zresm ! local real 894 CHARACTER(len=20) :: clname 895 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence 896 !!---------------------------------------------------------------------- 897 898 ! create file 899 IF( kt == nit000 .AND. kiter == 1 ) THEN 900 ! 901 IF( lwp ) THEN 902 WRITE(numout,*) 903 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 904 WRITE(numout,*) '~~~~~~~' 905 ENDIF 906 ! 907 IF( lwm ) THEN 908 clname = 'ice_cvg.nc' 909 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 910 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 911 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 912 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 913 istatus = NF90_ENDDEF(ncvgid) 914 ENDIF 915 ! 916 ENDIF 917 918 ! time 919 it = ( kt - 1 ) * kitermax + kiter 920 921 ! convergence 922 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 923 zresm = 0._wp 924 ELSE 925 DO jj = 1, jpj 926 DO ji = 1, jpi 927 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 928 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 929 END DO 930 END DO 931 zresm = MAXVAL( zres ) 932 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 933 ENDIF 934 935 IF( lwm ) THEN 936 ! write variables 937 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 938 ! close file 939 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid) 940 ENDIF 941 942 END SUBROUTINE rhg_cvg 855 943 856 944 … … 910 998 END SUBROUTINE rhg_evp_rst 911 999 1000 912 1001 #else 913 1002 !!---------------------------------------------------------------------- -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/iceistate.F90
r11715 r12937 41 41 ! !! ** namelist (namini) ** 42 42 LOGICAL, PUBLIC :: ln_iceini !: Ice initialization or not 43 LOGICAL, PUBLIC :: ln_iceini_file !: Ice initialization from 2D netcdf file 43 INTEGER, PUBLIC :: nn_iceini_file !: Ice initialization: 44 ! 0 = Initialise sea ice based on SSTs 45 ! 1 = Initialise sea ice from single category netcdf file 46 ! 2 = Initialise sea ice from multi category restart file 44 47 REAL(wp) :: rn_thres_sst 45 48 REAL(wp) :: rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n 46 49 REAL(wp) :: rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s 47 REAL(wp) :: rn_apd_ini_n, rn_hpd_ini_n 48 REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s 50 REAL(wp) :: rn_apd_ini_n, rn_hpd_ini_n, rn_hld_ini_n 51 REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s, rn_hld_ini_s 49 52 ! 50 ! ! if ln_iceini_file = T51 INTEGER , PARAMETER :: jpfldi = 9! maximum number of files to read53 ! ! if nn_iceini_file = 1 54 INTEGER , PARAMETER :: jpfldi = 10 ! maximum number of files to read 52 55 INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m) 53 56 INTEGER , PARAMETER :: jp_hts = 2 ! index of snw thickness (m) … … 59 62 INTEGER , PARAMETER :: jp_apd = 8 ! index of pnd fraction (-) 60 63 INTEGER , PARAMETER :: jp_hpd = 9 ! index of pnd depth (m) 64 INTEGER , PARAMETER :: jp_hld = 10 ! index of pnd lid depth (m) 61 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read) 62 66 ! … … 98 102 REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file 99 103 REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 100 REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini 104 REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini, zhlid_ini !data from namelist or nc file 101 105 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays 102 106 !! 103 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d 107 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d 104 108 !-------------------------------------------------------------------- 105 109 … … 155 159 a_ip (:,:,:) = 0._wp 156 160 v_ip (:,:,:) = 0._wp 157 a_ip_frac(:,:,:) = 0._wp 161 v_il (:,:,:) = 0._wp 162 a_ip_eff (:,:,:) = 0._wp 158 163 h_ip (:,:,:) = 0._wp 164 h_il (:,:,:) = 0._wp 159 165 ! 160 166 ! ice velocities … … 167 173 IF( ln_iceini ) THEN 168 174 ! !---------------! 169 IF( ln_iceini_file )THEN! Read a file !175 IF( nn_iceini_file == 1 )THEN ! Read a file ! 170 176 ! !---------------! 171 177 WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp … … 218 224 IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 219 225 & si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 226 ! 227 ! pond lid depth 228 IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) & 229 & si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 230 ! 220 231 zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) 232 zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) 221 233 ! 222 234 ! change the switch for the following … … 243 255 zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 244 256 zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 257 zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:) 245 258 ELSEWHERE 246 259 zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) … … 253 266 zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 254 267 zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 268 zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:) 255 269 END WHERE 256 270 ! … … 261 275 zapnd_ini(:,:) = 0._wp 262 276 zhpnd_ini(:,:) = 0._wp 277 zhlid_ini(:,:) = 0._wp 278 ENDIF 279 280 IF ( .NOT.ln_pnd_lids ) THEN 281 zhlid_ini(:,:) = 0._wp 263 282 ENDIF 264 283 … … 287 306 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini ) 288 307 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini ) 308 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti) , zhlid_ini ) 289 309 290 310 ! allocate temporary arrays 291 ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 292 & zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 311 ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), & 312 & zti_2d (npti,jpl), zts_2d (npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), & 313 & zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_2d(npti,jpl) ) 293 314 294 315 ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 295 CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), & 296 & zhi_2d , zhs_2d , zai_2d , & 297 & t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), & 298 & zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d ) 316 CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), & 317 & zhi_2d , zhs_2d , zai_2d , & 318 & t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), & 319 & s_i_1d(1:npti) , a_ip_1d(1:npti) , h_ip_1d(1:npti), h_il_1d(1:npti), & 320 & zti_2d , zts_2d , ztsu_2d , & 321 & zsi_2d , zaip_2d , zhip_2d , zhil_2d ) 299 322 300 323 ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) … … 312 335 CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip ) 313 336 CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip ) 337 CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d , h_il ) 314 338 315 339 ! deallocate temporary arrays 316 340 DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 317 & zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )341 & zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d ) 318 342 319 343 ! calculate extensive and intensive variables … … 357 381 358 382 ! Melt ponds 359 WHERE( a_i > epsi10 ) 360 a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 361 ELSEWHERE 362 a_ip_frac(:,:,:) = 0._wp 383 WHERE( a_i > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 384 ELSEWHERE ; a_ip_eff(:,:,:) = 0._wp 363 385 END WHERE 364 386 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 387 v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:) 365 388 366 389 ! specific temperatures for coupled runs … … 434 457 e_s_b (:,:,:,:) = e_s (:,:,:,:) 435 458 sv_i_b (:,:,:) = sv_i (:,:,:) 436 oa_i_b (:,:,:) = oa_i (:,:,:)437 459 u_ice_b(:,:) = u_ice(:,:) 438 460 v_ice_b(:,:) = v_ice(:,:) … … 463 485 ! 464 486 CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files 465 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd 487 TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd, sn_hld 466 488 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 467 489 ! 468 NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &490 NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, & 469 491 & rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, & 470 492 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, & 471 493 & rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, & 472 & rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &473 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir494 & rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, rn_hld_ini_n, rn_hld_ini_s, & 495 & sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, sn_hld, cn_dir 474 496 !!----------------------------------------------------------------------------- 475 497 ! … … 485 507 slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi 486 508 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms 487 slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd 509 slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd ; slf_i(jp_hld) = sn_hld 488 510 ! 489 511 IF(lwp) THEN ! control print … … 493 515 WRITE(numout,*) ' Namelist namini:' 494 516 WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini 495 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file517 WRITE(numout,*) ' ice initialization from a netcdf file nn_iceini_file = ', nn_iceini_file 496 518 WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 497 IF( ln_iceini .AND. .NOT.ln_iceini_file) THEN519 IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN 498 520 WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s 499 521 WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s … … 505 527 WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s 506 528 WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s 529 WRITE(numout,*) ' initial pnd lid depth in the north-south rn_hld_ini = ', rn_hld_ini_n,rn_hld_ini_s 507 530 ENDIF 508 531 ENDIF 509 532 ! 510 IF( ln_iceini_file) THEN ! Ice initialization using input file533 IF( nn_iceini_file == 1 ) THEN ! Ice initialization using input file 511 534 ! 512 535 ! set si structure … … 529 552 rn_apd_ini_n = 0. ; rn_apd_ini_s = 0. 530 553 rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0. 531 CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' ) 554 rn_hld_ini_n = 0. ; rn_hld_ini_s = 0. 555 CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 & rn_hld_ini = 0 when no ponds' ) 556 ENDIF 557 ! 558 IF( .NOT.ln_pnd_lids ) THEN 559 rn_hld_ini_n = 0. ; rn_hld_ini_s = 0. 532 560 ENDIF 533 561 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/iceitd.F90
r11715 r12937 411 411 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 412 412 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 413 CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il ) 413 414 CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 414 415 DO jl = 1, jpl … … 483 484 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 484 485 v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 486 ! 487 IF ( ln_pnd_lids ) THEN ! Pond lid volume 488 ztrans = v_il_2d(ji,jl1) * zworka(ji) 489 v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans 490 v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans 491 ENDIF 485 492 ENDIF 486 493 ! … … 527 534 ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 528 535 ! because of truncation error ( i.e. 1. - 1. /= 0 ) 529 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )536 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d ) 530 537 531 538 ! at_i must be <= rn_amax … … 555 562 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 556 563 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 564 CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il ) 557 565 CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 558 566 DO jl = 1, jpl -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icerst.F90
r11715 r12937 132 132 CALL iom_rstput( iter, nitrst, numriw, 'a_ip' , a_ip ) 133 133 CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip ) 134 CALL iom_rstput( iter, nitrst, numriw, 'v_il' , v_il ) 134 135 ! Snow enthalpy 135 136 DO jk = 1, nlay_s … … 171 172 INTEGER :: jk 172 173 LOGICAL :: llok 173 INTEGER :: id0, id1, id2, id3, id4 ! local integer174 INTEGER :: id0, id1, id2, id3, id4, id5 ! local integer 174 175 CHARACTER(len=25) :: znam 175 176 CHARACTER(len=2) :: zchar, zchar1 … … 250 251 v_ip(:,:,:) = 0._wp 251 252 ENDIF 253 ! melt pond lids 254 id3 = iom_varid( numrir, 'v_il' , ldstop = .FALSE. ) 255 IF( id3 > 0 ) THEN 256 CALL iom_get( numrir, jpdom_autoglo, 'v_il', v_il) 257 ELSE 258 IF(lwp) WRITE(numout,*) ' ==>> previous run without melt ponds lids output then set it to zero' 259 v_il(:,:,:) = 0._wp 260 ENDIF 252 261 ! fields needed for Met Office (Jules) coupling 253 262 IF( ln_cpl ) THEN 254 id 3= iom_varid( numrir, 'cnd_ice' , ldstop = .FALSE. )255 id 4= iom_varid( numrir, 't1_ice' , ldstop = .FALSE. )256 IF( id 3 > 0 .AND. id4> 0 ) THEN ! fields exist263 id4 = iom_varid( numrir, 'cnd_ice' , ldstop = .FALSE. ) 264 id5 = iom_varid( numrir, 't1_ice' , ldstop = .FALSE. ) 265 IF( id4 > 0 .AND. id5 > 0 ) THEN ! fields exist 257 266 CALL iom_get( numrir, jpdom_autoglo, 'cnd_ice', cnd_ice ) 258 267 CALL iom_get( numrir, jpdom_autoglo, 't1_ice' , t1_ice ) … … 274 283 CALL ice_istate( nit000 ) 275 284 ! 276 IF( .NOT.ln_iceini .OR. .NOT.ln_iceini_file) &277 & CALL ctl_stop('STOP', 'ice_rst_read: you need ln_ice_ini=T and ln_iceini_file=T')285 IF( .NOT.ln_iceini .OR. nn_iceini_file == 0 ) & 286 & CALL ctl_stop('STOP', 'ice_rst_read: you need ln_ice_ini=T and nn_iceini_file=0') 278 287 ! 279 288 ENDIF -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icesbc.F90
r11715 r12937 116 116 INTEGER :: ji, jj, jl ! dummy loop index 117 117 REAL(wp) :: zmiss_val ! missing value retrieved from xios 118 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb_os, zalb_cs ! ice albedo under overcast/clear sky 119 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zalb, zmsk00 ! 2D workspace 118 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zalb, zmsk00 ! 2D workspace 120 119 !!-------------------------------------------------------------------- 121 120 ! … … 131 130 CALL iom_miss_val( "icetemp", zmiss_val ) 132 131 133 ! --- cloud-sky and overcast-sky ice albedos --- ! 134 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) 135 136 ! albedo depends on cloud fraction because of non-linear spectral effects 137 !!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument ! 138 alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 139 ! 132 ! --- ice albedo --- ! 133 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, cloud_fra, alb_ice ) 134 140 135 ! 141 136 SELECT CASE( ksbc ) !== fluxes over sea ice ==! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icestp.F90
r11715 r12937 252 252 ! 253 253 ! ! Initial sea-ice state 254 IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst 254 255 IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN 256 CALL ice_rst_read ! start from a restart file 257 ELSE 255 258 CALL ice_istate_init 256 CALL ice_istate( nit000 ) 257 ELSE ! start from a restart file 258 CALL ice_rst_read 259 CALL ice_istate( nit000 ) ! start from rest: sea-ice deduced from sst 259 260 ENDIF 260 261 CALL ice_var_glo2eqv … … 363 364 v_s_b (:,:,:) = v_s (:,:,:) ! snow volume 364 365 sv_i_b(:,:,:) = sv_i(:,:,:) ! salt content 365 oa_i_b(:,:,:) = oa_i(:,:,:) ! areal age content366 366 e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy 367 367 e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy … … 372 372 h_i_b(:,:,:) = 0._wp 373 373 h_s_b(:,:,:) = 0._wp 374 END WHERE375 376 WHERE( a_ip(:,:,:) >= epsi20 )377 h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) ! ice pond thickness378 ELSEWHERE379 h_ip_b(:,:,:) = 0._wp380 374 END WHERE 381 375 ! … … 421 415 hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp 422 416 hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp 423 hfx_err_rem(:,:) = 0._wp424 417 hfx_err_dif(:,:) = 0._wp 425 418 wfx_err_sub(:,:) = 0._wp -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icethd.F90
r11715 r12937 51 51 LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F) 52 52 LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F) 53 LOGICAL :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean 53 54 54 55 !! * Substitutions … … 164 165 ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 165 166 IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 166 fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 167 IF( ln_leadhfx ) THEN ; fhld(ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 168 ELSE ; fhld(ji,jj) = 0._wp 169 ENDIF 167 170 qlead(ji,jj) = 0._wp 168 171 ELSE … … 354 357 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) ) 355 358 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 356 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )359 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) ) 357 360 ! 358 361 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice ) … … 406 409 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 407 410 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 408 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )409 411 CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai ) 410 412 ! … … 441 443 sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 442 444 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 445 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 443 446 oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti) 444 447 … … 460 463 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) ) 461 464 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 462 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )465 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) ) 463 466 ! 464 467 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) … … 498 501 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 499 502 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 500 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )501 503 CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai ) 502 504 ! … … 515 517 CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 516 518 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) ) 519 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) ) 517 520 CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) ) 518 521 ! … … 536 539 INTEGER :: ios ! Local integer output status for namelist read 537 540 !! 538 NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS 541 NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS, ln_leadhfx 539 542 !!------------------------------------------------------------------- 540 543 ! … … 552 555 WRITE(numout,*) '~~~~~~~~~~~~' 553 556 WRITE(numout,*) ' Namelist namthd:' 554 WRITE(numout,*) ' activate ice thick change from top/bot (T) or not (F) ln_icedH = ', ln_icedH 555 WRITE(numout,*) ' activate lateral melting (T) or not (F) ln_icedA = ', ln_icedA 556 WRITE(numout,*) ' activate ice growth in open-water (T) or not (F) ln_icedO = ', ln_icedO 557 WRITE(numout,*) ' activate gravity drainage and flushing (T) or not (F) ln_icedS = ', ln_icedS 557 WRITE(numout,*) ' activate ice thick change from top/bot (T) or not (F) ln_icedH = ', ln_icedH 558 WRITE(numout,*) ' activate lateral melting (T) or not (F) ln_icedA = ', ln_icedA 559 WRITE(numout,*) ' activate ice growth in open-water (T) or not (F) ln_icedO = ', ln_icedO 560 WRITE(numout,*) ' activate gravity drainage and flushing (T) or not (F) ln_icedS = ', ln_icedS 561 WRITE(numout,*) ' heat in the leads is used to melt sea-ice before warming the ocean ln_leadhfx = ', ln_leadhfx 558 562 ENDIF 559 563 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icethd_ent.F90
r11715 r12937 128 128 ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_do), 129 129 ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 130 DO ji = 1, npti131 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice * &132 & ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) )133 END DO130 !DO ji = 1, npti 131 ! hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice * & 132 ! & ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 133 !END DO 134 134 135 135 END SUBROUTINE ice_thd_ent -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icethd_pnd.F90
r11715 r12937 88 88 ! 89 89 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 90 a_ip_frac_1d(ji) = rn_apnd91 90 h_ip_1d(ji) = rn_hpnd 92 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 91 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 92 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 93 93 ELSE 94 a_ip_frac_1d(ji) = 0._wp95 94 h_ip_1d(ji) = 0._wp 96 95 a_ip_1d(ji) = 0._wp 96 h_il_1d(ji) = 0._wp 97 97 ENDIF 98 98 ! … … 106 106 !! *** ROUTINE pnd_H12 *** 107 107 !! 108 !! ** Purpose : Compute melt pond evolution 109 !! 110 !! ** Method : Empirical method. A fraction of meltwater is accumulated in ponds 111 !! and sent to ocean when surface is freezing 112 !! 113 !! pond growth: Vp = Vp + dVmelt 114 !! with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i 115 !! pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) 116 !! with Tp = -2degC 117 !! 118 !! ** Tunable parameters : (no real expertise yet, ideas?) 108 !! ** Purpose : Compute melt pond evolution 109 !! 110 !! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing 111 !! We work with volumes and then redistribute changes into thickness and concentration 112 !! assuming linear relationship between the two. 113 !! 114 !! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 --- 115 !! dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i 116 !! dh_i = meltwater from ice surface melt 117 !! dh_s = meltwater from snow melt 118 !! (1-r) = fraction of melt water that is not flushed 119 !! 120 !! - limtations: a_ip must not exceed (1-r)*a_i 121 !! h_ip must not exceed 0.5*h_i 122 !! 123 !! - pond shrinking: 124 !! if lids: Vp = Vp -dH * a_ip 125 !! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 --- 126 !! 127 !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 128 !! H = lid thickness 129 !! Lf = latent heat of fusion 130 !! Tp = -2C 131 !! 132 !! And solved implicitely as: 133 !! H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0 134 !! 135 !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- 136 !! 137 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi --- from Flocco et al 2007 --- 138 !! perm = permability of sea-ice 139 !! visc = water viscosity 140 !! Hp = height of top of the pond above sea-level 141 !! Hi = ice thickness thru which there is flushing 142 !! 143 !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness 144 !! 145 !! - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip: 146 !! a_ip/a_i = a_ip_frac = h_ip / zaspect 147 !! 148 !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 119 149 !! 120 !! ** Note : Stolen from CICE for quick test of the melt pond 121 !! radiation and freshwater interfaces 122 !! Coupling can be radiative AND freshwater 123 !! Advection, ridging, rafting are called 124 !! 125 !! ** References : Holland, M. M. et al (J Clim 2012) 126 !!------------------------------------------------------------------- 127 REAL(wp), PARAMETER :: zrmin = 0.15_wp ! minimum fraction of available meltwater retained for melt ponding 128 REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum - - - - - 129 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 130 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 131 ! 132 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 133 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding 134 REAL(wp) :: z1_Tp ! inverse reference temperature 135 REAL(wp) :: z1_rhow ! inverse freshwater density 136 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 137 REAL(wp) :: zfac, zdum 138 ! 139 INTEGER :: ji ! loop indices 140 !!------------------------------------------------------------------- 141 z1_rhow = 1._wp / rhow 142 z1_zpnd_aspect = 1._wp / zpnd_aspect 143 z1_Tp = 1._wp / zTp 150 !! ** Note : mostly stolen from CICE 151 !! 152 !! ** References : Flocco and Feltham (JGR, 2007) 153 !! Flocco et al (JGR, 2010) 154 !! Holland et al (J. Clim, 2012) 155 !!------------------------------------------------------------------- 156 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array 157 !! 158 REAL(wp), PARAMETER :: zaspect = 0.8_wp ! pond aspect ratio 159 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 160 REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity 161 !! 162 REAL(wp) :: zfr_mlt, zdv_mlt ! fraction and volume of available meltwater retained for melt ponding 163 REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes 164 REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh 165 REAL(wp) :: zv_ip_max ! max pond volume allowed 166 REAL(wp) :: zdT ! zTp-t_su 167 REAL(wp) :: zsbr ! Brine salinity 168 REAL(wp) :: zperm ! permeability of sea ice 169 REAL(wp) :: zfac, zdum ! temporary arrays 170 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 171 !! 172 INTEGER :: ji, jk ! loop indices 173 !!------------------------------------------------------------------- 174 z1_rhow = 1._wp / rhow 175 z1_aspect = 1._wp / zaspect 176 z1_Tp = 1._wp / zTp 144 177 145 178 DO ji = 1, npti 146 ! !--------------------------------!147 IF( h_i_1d(ji) < rn_himin ) THEN ! Case ice thickness < rn_himin!148 ! !--------------------------------!149 !--- Remove ponds on thin ice 179 ! !----------------------------------------------------! 180 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 181 ! !----------------------------------------------------! 182 !--- Remove ponds on thin ice or tiny ice fractions 150 183 a_ip_1d(ji) = 0._wp 151 a_ip_frac_1d(ji) = 0._wp152 184 h_ip_1d(ji) = 0._wp 153 ! !--------------------------------! 154 ELSE ! Case ice thickness >= rn_himin ! 155 ! !--------------------------------! 156 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! record pond volume at previous time step 157 ! 158 ! available meltwater for melt ponding [m, >0] and fraction 159 zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 160 zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) ! from CICE doc 161 !zfr_mlt = zrmin + zrmax * a_i_1d(ji) ! from Holland paper 162 ! 163 !--- Pond gowth ---! 164 ! v_ip should never be negative, otherwise code crashes 165 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 166 ! 167 ! melt pond mass flux (<0) 185 h_il_1d(ji) = 0._wp 186 ! 187 ! clem: problem with conservation or not ? 188 ! !--------------------------------! 189 ELSE ! Case ice thickness >= rn_himin ! 190 ! !--------------------------------! 191 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 192 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 193 ! 194 !------------------! 195 ! case ice melting ! 196 !------------------! 197 ! 198 !--- available meltwater for melt ponding ---! 199 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 200 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) in H12 = fraction of melt water that is not flushed 201 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 202 ! 203 !--- overflow ---! 204 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 205 ! a_ip_max = zfr_mlt * a_i 206 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 207 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 208 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 209 210 ! If pond depth exceeds half the ice thickness then reduce the pond volume 211 ! h_ip_max = 0.5 * h_i 212 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 213 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 214 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 215 216 !--- Pond growing ---! 217 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 218 ! 219 !--- Lid melting ---! 220 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 221 ! 222 !--- mass flux ---! 168 223 IF( zdv_mlt > 0._wp ) THEN 169 zfac = z fr_mlt * zdv_mlt * rhow * r1_rdtice224 zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1] 170 225 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 171 226 ! 172 ! adjust ice/snow melting flux to balance melt pond flux (>0) 173 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 227 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 174 228 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 175 229 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 176 230 ENDIF 231 232 !-------------------! 233 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 234 !-------------------! 235 ! 236 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 177 237 ! 178 238 !--- Pond contraction (due to refreezing) ---! 179 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 180 ! 181 ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 182 ! h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i 183 a_ip_1d(ji) = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) ) 184 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 185 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 239 IF( ln_pnd_lids ) THEN 240 ! 241 !--- Lid growing and subsequent pond shrinking ---! 242 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 243 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 244 245 ! Lid growing 246 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 247 248 ! Pond shrinking 249 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 250 251 ELSE 252 ! Pond shrinking 253 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 254 ENDIF 255 ! 256 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 257 ! v_ip = h_ip * a_ip 258 ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 259 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 260 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 261 262 !---------------! 263 ! Pond flushing ! 264 !---------------! 265 IF( ln_pnd_flush ) THEN 266 ! height of top of the pond above sea-level 267 zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 268 269 ! Calculate the permeability of the ice (Assur 1958) 270 DO jk = 1, nlay_i 271 zsbr = - 1.2_wp & 272 & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 273 & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 274 & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 ! clem: error here the factor was 0.01878 instead of 0.0178 (cf Flocco 2010) 275 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 276 END DO 277 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 278 279 ! Do the drainage using Darcy's law 280 zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 281 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 282 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 283 284 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 285 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 286 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 287 288 ENDIF 289 290 !--- Corrections and lid thickness ---! 291 IF( ln_pnd_lids ) THEN 292 !--- retrieve lid thickness from volume ---! 293 IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 294 ELSE ; h_il_1d(ji) = 0._wp 295 ENDIF 296 !--- remove ponds if lids are much larger than ponds ---! 297 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 298 a_ip_1d(ji) = 0._wp 299 h_ip_1d(ji) = 0._wp 300 h_il_1d(ji) = 0._wp 301 ENDIF 302 ENDIF 186 303 ! 187 304 ENDIF 305 188 306 END DO 189 307 ! … … 205 323 INTEGER :: ios, ioptio ! Local integer 206 324 !! 207 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 325 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_lids, ln_pnd_flush, rn_apnd_min, rn_apnd_max, & 326 & ln_pnd_CST, rn_apnd, rn_hpnd, & 327 & ln_pnd_alb 208 328 !!------------------------------------------------------------------- 209 329 ! … … 221 341 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 222 342 WRITE(numout,*) ' Namelist namicethd_pnd:' 223 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 224 WRITE(numout,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 225 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST 226 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd 227 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd 228 WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb 343 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 344 WRITE(numout,*) ' Evolutive melt pond fraction and depth ln_pnd_H12 = ', ln_pnd_H12 345 WRITE(numout,*) ' Melt ponds can have frozen lids ln_pnd_lids = ', ln_pnd_lids 346 WRITE(numout,*) ' Allow ponds to flush thru the ice ln_pnd_flush = ', ln_pnd_flush 347 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min 348 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max 349 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST 350 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd 351 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd 352 WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb 229 353 ENDIF 230 354 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/iceupdate.F90
r11715 r12937 94 94 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 95 95 REAL(wp) :: zqsr ! New solar flux received by the ocean 96 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 97 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb_cs, zalb_os ! 3D workspace 96 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 98 97 !!--------------------------------------------------------------------- 99 98 IF( ln_timing ) CALL timing_start('ice_update') … … 185 184 ! Snow/ice albedo (only if sent to coupler, useless in forced mode) 186 185 !------------------------------------------------------------------ 187 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 188 ! 189 alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 186 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, cloud_fra, alb_ice ) ! ice albedo 187 190 188 ! 191 189 IF( lrst_ice ) THEN !* write snwice_mass fields in the restart file -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icevar.F90
r11715 r12937 113 113 at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds 114 114 vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) 115 vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) 115 116 ! 116 117 ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction … … 161 162 ! 162 163 ! ! mean melt pond depth 163 WHERE( at_ip(:,:) > epsi20 ) ; hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) 164 ELSEWHERE ; hm_ip(:,:) = 0._wp 164 WHERE( at_ip(:,:) > epsi20 ) ; hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) ; hm_il(:,:) = vt_il(:,:) / at_ip(:,:) 165 ELSEWHERE ; hm_ip(:,:) = 0._wp ; hm_il(:,:) = 0._wp 165 166 END WHERE 166 167 ! … … 184 185 REAL(wp) :: zhmax, z1_zhmax ! - - 185 186 REAL(wp) :: zlay_i, zlay_s ! - - 186 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i 187 REAL(wp), PARAMETER :: zhl_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 188 REAL(wp), PARAMETER :: zhl_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation 189 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i, z1_a_ip 187 190 !!------------------------------------------------------------------- 188 191 … … 203 206 ELSEWHERE ; z1_v_i(:,:,:) = 0._wp 204 207 END WHERE 208 ! 209 WHERE( a_ip(:,:,:) > epsi20 ) ; z1_a_ip(:,:,:) = 1._wp / a_ip(:,:,:) 210 ELSEWHERE ; z1_a_ip(:,:,:) = 0._wp 211 END WHERE 205 212 ! !--- ice thickness 206 213 h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) … … 217 224 ! !--- ice age 218 225 o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) 219 ! !--- pond fraction and thickness 226 ! !--- pond and lid thickness 227 h_ip(:,:,:) = v_ip(:,:,:) * z1_a_ip(:,:,:) 228 h_il(:,:,:) = v_il(:,:,:) * z1_a_ip(:,:,:) 229 ! !--- melt pond effective area (used for albedo) 220 230 a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:) 221 WHERE( a_ip_frac(:,:,:) > epsi20 ) ; h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:) 222 ELSEWHERE ; h_ip(:,:,:) = 0._wp 223 END WHERE 224 ! 231 WHERE ( h_il(:,:,:) <= zhl_min ) ; a_ip_eff(:,:,:) = a_ip_frac(:,:,:) ! lid is very thin. Expose all the pond 232 ELSEWHERE( h_il(:,:,:) >= zhl_max ) ; a_ip_eff(:,:,:) = 0._wp ! lid is very thick. Cover all the pond up with ice and snow 233 ELSEWHERE ; a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * & ! lid is in between. Expose part of the pond 234 & ( h_il(:,:,:) - zhl_min ) / ( zhl_max - zhl_min ) 235 END WHERE 225 236 ! !--- salinity (with a minimum value imposed everywhere) 226 237 IF( nn_icesal == 2 ) THEN … … 289 300 sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:) 290 301 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 302 v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:) 291 303 ! 292 304 END SUBROUTINE ice_var_eqv2glo … … 533 545 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 534 546 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 547 v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) 535 548 ! 536 549 END DO … … 555 568 556 569 557 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )570 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 558 571 !!------------------------------------------------------------------- 559 572 !! *** ROUTINE ice_var_zapneg *** … … 570 583 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 571 584 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 585 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid volume 572 586 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 573 587 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 636 650 WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp 637 651 WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 638 !but it does not change conservation, so keep it this way is ok652 WHERE( pv_il (:,:,:) < 0._wp ) pv_il (:,:,:) = 0._wp ! but it does not change conservation, so keep it this way is ok 639 653 ! 640 654 END SUBROUTINE ice_var_zapneg 641 655 642 656 643 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, p e_s, pe_i )657 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 644 658 !!------------------------------------------------------------------- 645 659 !! *** ROUTINE ice_var_roundoff *** … … 654 668 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 655 669 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume 670 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_il ! melt pond lid volume 656 671 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content 657 672 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content … … 668 683 WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 669 684 WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 685 IF( ln_pnd_lids ) THEN 686 WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:) = 0._wp ! v_il must be >= 0 687 ENDIF 670 688 ENDIF 671 689 ! … … 786 804 !! ** Purpose : converting N-cat ice to jpl ice categories 787 805 !!------------------------------------------------------------------- 788 SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, &789 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip)806 SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, & 807 & ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 790 808 !!------------------------------------------------------------------- 791 809 !! ** Purpose : converting 1-cat ice to 1 ice category … … 793 811 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 794 812 REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 795 REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds796 REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds813 REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds 814 REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ! output ice/snow temp & sal & ponds 797 815 !!------------------------------------------------------------------- 798 816 ! == thickness and concentration == ! … … 808 826 pa_ip(:) = patip(:) 809 827 ph_ip(:) = phtip(:) 828 ph_il(:) = phtil(:) 810 829 811 830 END SUBROUTINE ice_var_itd_1c1c 812 831 813 SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, &814 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip)832 SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, & 833 & ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 815 834 !!------------------------------------------------------------------- 816 835 !! ** Purpose : converting N-cat ice to 1 ice category … … 818 837 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 819 838 REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 820 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds821 REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds839 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds 840 REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ! output ice/snow temp & sal & ponds 822 841 ! 823 842 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs … … 854 873 ! == ponds == ! 855 874 pa_ip(:) = SUM( patip(:,:), dim=2 ) 856 WHERE( pa_ip(:) /= 0._wp ) ; ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 857 ELSEWHERE ; ph_ip(:) = 0._wp 875 WHERE( pa_ip(:) /= 0._wp ) 876 ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 877 ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 878 ELSEWHERE 879 ph_ip(:) = 0._wp 880 ph_il(:) = 0._wp 858 881 END WHERE 859 882 ! … … 862 885 END SUBROUTINE ice_var_itd_Nc1c 863 886 864 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, &865 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip)887 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, & 888 & ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 866 889 !!------------------------------------------------------------------- 867 890 !! … … 885 908 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 886 909 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 887 REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds888 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds910 REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds 911 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ! output ice/snow temp & sal & ponds 889 912 ! 890 913 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti … … 997 1020 END WHERE 998 1021 END DO 1022 ! keep the same v_il/v_i ratio for each category 1023 WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) ) 1024 ELSEWHERE ; zfra(:) = 0._wp 1025 END WHERE 1026 DO jl = 1, jpl 1027 WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 1028 ELSEWHERE ; ph_il(:,jl) = 0._wp 1029 END WHERE 1030 END DO 999 1031 DEALLOCATE( zfra ) 1000 1032 ! 1001 1033 END SUBROUTINE ice_var_itd_1cMc 1002 1034 1003 SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, &1004 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip)1035 SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, & 1036 & ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) 1005 1037 !!------------------------------------------------------------------- 1006 1038 !! … … 1033 1065 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 1034 1066 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 1035 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds1036 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds1067 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip, phtil ! input ice/snow temp & sal & ponds 1068 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ! output ice/snow temp & sal & ponds 1037 1069 ! 1038 1070 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2 … … 1063 1095 pa_ip(:,:) = patip(:,:) 1064 1096 ph_ip(:,:) = phtip(:,:) 1097 ph_il(:,:) = phtil(:,:) 1065 1098 ! ! ---------------------- ! 1066 1099 ELSEIF( icat == 1 ) THEN ! input cat = 1 ! … … 1068 1101 CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 1069 1102 & ph_i(:,:), ph_s(:,:), pa_i (:,:), & 1070 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), &1071 & pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:) )1103 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), & 1104 & pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:) ) 1072 1105 ! ! ---------------------- ! 1073 1106 ELSEIF( jpl == 1 ) THEN ! output cat = 1 ! … … 1075 1108 CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 1076 1109 & ph_i(:,1), ph_s(:,1), pa_i (:,1), & 1077 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), &1078 & pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1) )1110 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), & 1111 & pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1) ) 1079 1112 ! ! ----------------------- ! 1080 1113 ELSE ! input cat /= output cat ! … … 1218 1251 END WHERE 1219 1252 END DO 1253 ! keep the same v_il/v_i ratio for each category 1254 WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 1255 zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 1256 ELSEWHERE 1257 zfra(:) = 0._wp 1258 END WHERE 1259 DO jl = 1, jpl 1260 WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 1261 ELSEWHERE ; ph_il(:,jl) = 0._wp 1262 END WHERE 1263 END DO 1220 1264 DEALLOCATE( zfra ) 1221 1265 ! -
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icewri.F90
r11715 r12937 116 116 IF( iom_use('icehpnd' ) ) CALL iom_put( 'icehpnd', hm_ip * zmsk00 ) ! melt pond depth 117 117 IF( iom_use('icevpnd' ) ) CALL iom_put( 'icevpnd', vt_ip * zmsk00 ) ! melt pond total volume per unit area 118 IF( iom_use('icehlid' ) ) CALL iom_put( 'icehlid', hm_il * zmsk00 ) ! melt pond lid depth 119 IF( iom_use('icevlid' ) ) CALL iom_put( 'icevlid', vt_il * zmsk00 ) ! melt pond lid total volume per unit area 118 120 ! salt 119 121 IF( iom_use('icesalt' ) ) CALL iom_put( 'icesalt', sm_i * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity … … 149 151 150 152 ! --- category-dependent fields --- ! 153 IF( iom_use('icehlid_cat' ) ) CALL iom_put( 'icehlid_cat' , h_il * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories 154 IF( iom_use('iceaepnd_cat') ) CALL iom_put( 'iceaepnd_cat', a_ip_eff * zmsk00l ) ! melt pond effective frac for categories 151 155 IF( iom_use('icemask_cat' ) ) CALL iom_put( 'icemask_cat' , zmsk00l ) ! ice mask 0% 152 156 IF( iom_use('iceconc_cat' ) ) CALL iom_put( 'iceconc_cat' , a_i * zmsk00l ) ! area for categories … … 162 166 IF( iom_use('icebrv_cat' ) ) CALL iom_put( 'icebrv_cat' , bv_i * 100. * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume 163 167 IF( iom_use('iceapnd_cat' ) ) CALL iom_put( 'iceapnd_cat' , a_ip * zmsk00l ) ! melt pond frac for categories 164 IF( iom_use('icehpnd_cat' ) ) CALL iom_put( 'icehpnd_cat' , h_ip * zmsk00l + zmiss_val * ( 1._wp - zmsk00l )) ! melt pond frac for categories168 IF( iom_use('icehpnd_cat' ) ) CALL iom_put( 'icehpnd_cat' , h_ip * zmsk00l ) ! melt pond frac for categories 165 169 IF( iom_use('iceafpnd_cat') ) CALL iom_put( 'iceafpnd_cat', a_ip_frac * zmsk00l ) ! melt pond frac for categories 166 IF( iom_use('icealb_cat' ) ) CALL iom_put( 'icealb_cat' , alb_ice * zmsk00l + zmiss_val * ( 1._wp - zmsk00l )) ! ice albedo for categories170 IF( iom_use('icealb_cat' ) ) CALL iom_put( 'icealb_cat' , alb_ice * zmsk00l ) ! ice albedo for categories 167 171 168 172 !------------------
Note: See TracChangeset
for help on using the changeset viewer.