Changeset 8623
- Timestamp:
- 2017-10-12T23:24:57+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 7 edited
-
ice1D.F90 (modified) (4 diffs)
-
icectl.F90 (modified) (3 diffs)
-
icethd.F90 (modified) (7 diffs)
-
icethd_dh.F90 (modified) (5 diffs)
-
icethd_pnd.F90 (modified) (6 diffs)
-
iceupdate.F90 (modified) (1 diff)
-
icevar.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
r8565 r8623 89 89 90 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sprecip_1d !: <==> the 2D sprecip 91 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: at_i_1d 91 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: at_i_1d !: <==> the 2D at_i 92 92 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ato_i_1d !: <==> the 2D ato_i 93 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhtur_1d !: <==> the 2D fhtur … … 103 103 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_i_1d !: <==> the 2D a_i 104 104 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ib_1d !: <==> the 2D a_i_b 105 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_i_1d !:106 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_ib_1d !:107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_s_1d !:105 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_i_1d !: 106 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_ib_1d !: 107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_s_1d !: 108 108 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_su !: Surface Conduction flux 109 109 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_bo_i !: Bottom Conduction flux … … 112 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_sub !: Ice surface sublimation [m] 113 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_bott !: Ice bottom accretion/ablation [m] 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_s_mlt !: Snow melt [m] 114 115 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_snowice !: Snow ice formation [m of ice] 115 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_1d !: Ice bulk salinity [ppt]116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_1d !: Ice bulk salinity [ppt] 116 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_new !: Salinity of new ice at the bottom 118 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_i_1d !: 119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_s_1d !: 120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sv_i_1d !: 121 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: oa_i_1d !: 117 122 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_1d !: 118 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_ip_1d !: 119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_i_1d !: 120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_s_1d !: 121 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sv_i_1d !: 122 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: oa_i_1d !: 123 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_i_1d !: corresponding to the 2D var t_i 126 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sz_i_1d !: profiled ice salinity 127 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e_i_1d !: Ice enthalpy per unit volume 128 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e_s_1d !: Snow enthalpy per unit volume 129 130 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: eh_i_old !: ice heat content (q*h, J.m-2) 131 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_old !: ice thickness layer (m) 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: h_ip_1d !: 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_frac_1d !: 126 127 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s 128 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_i_1d !: corresponding to the 2D var t_i 129 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sz_i_1d !: profiled ice salinity 130 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e_i_1d !: Ice enthalpy per unit volume 131 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e_s_1d !: Snow enthalpy per unit volume 132 133 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: eh_i_old !: ice heat content (q*h, J.m-2) 134 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h_i_old !: ice thickness layer (m) 132 135 133 136 ! Conduction flux diagnostics (SIMIP) … … 200 203 & h_i_1d (jpij) , h_ib_1d (jpij) , h_s_1d (jpij) , fc_su (jpij) , fc_bo_i(jpij) , & 201 204 & dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) , & 202 & dh_i_bott(jpij) , dh_s nowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , &205 & dh_i_bott(jpij) , dh_s_mlt(jpij), dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , & 203 206 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , & 207 & h_ip_1d (jpij) , a_ip_frac_1d(jpij) , & 204 208 & sv_i_1d (jpij) , oa_i_1d (jpij) , STAT=ierr(ii) ) 205 209 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
r8567 r8623 76 76 IF( icount == 0 ) THEN 77 77 ! ! water flux 78 pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 79 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_ice_sub(:,:) + & 80 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:) & 78 pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 79 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:) + & 80 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + & 81 & wfx_ice_sub(:,:) + wfx_spr(:,:) & 81 82 & ) * e1e2t(:,:) ) * zconv 82 83 ! … … 96 97 97 98 pdiag_t = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 98 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 99 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv 99 100 100 101 ELSEIF( icount == 1 ) THEN 101 102 102 103 ! water flux 103 zfv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 104 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_ice_sub(:,:) + & 105 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:) & 104 zfv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 105 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:) + & 106 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + & 107 & wfx_ice_sub(:,:) + wfx_spr(:,:) & 106 108 & ) * e1e2t(:,:) ) * zconv - pdiag_fv 107 109 … … 117 119 118 120 ! outputs 119 zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t 121 zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv & 120 122 & - pdiag_v ) * r1_rdtice - zfv ) * rday 121 123 122 zs = ( ( glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t ) * zconv &124 zs = ( ( glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t ) * zconv & 123 125 & - pdiag_s ) * r1_rdtice + zfs ) * rday 124 126 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
r8598 r8623 231 231 s_i_new (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp ! --- some init --- ! (important to have them here) 232 232 dh_i_surf (1:npti) = 0._wp ; dh_i_bott(1:npti) = 0._wp 233 dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp 233 dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 234 234 ! 235 235 IF( ln_icedH ) THEN ! --- growing/melting --- ! 236 236 CALL ice_thd_zdf ! Ice/Snow Temperature profile 237 237 CALL ice_thd_dh ! Ice/Snow thickness 238 CALL ice_thd_pnd ! Melt ponds formation 238 239 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 239 240 ENDIF … … 273 274 ! 274 275 IF( ln_icedO ) CALL ice_thd_do ! --- frazil ice growing in leads --- ! 275 !276 CALL ice_thd_pnd( kt ) ! --- Melt ponds --- !277 276 ! 278 277 ! controls … … 365 364 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i ) 366 365 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl) ) 367 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,kl) )368 CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s(:,:,kl) )366 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl) ) 367 CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl) ) 369 368 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl) ) 370 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i(:,:,kl) )369 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl) ) 371 370 DO jk = 1, nlay_s 372 CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )373 CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )371 CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) ) 372 CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) ) 374 373 END DO 375 374 DO jk = 1, nlay_i 376 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) ) 377 CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) ) 378 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 379 END DO 375 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl) ) 376 CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl) ) 377 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 378 END DO 379 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) ) 380 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 381 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 380 382 ! 381 383 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice ) … … 406 408 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 407 409 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 410 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 408 411 ! 409 412 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) … … 457 460 ! 458 461 ! Change thickness to volume (replaces routine ice_var_eqv2glo) 459 v_i_1d(1:npti) = h_i_1d(1:npti) * a_i_1d(1:npti) 460 v_s_1d(1:npti) = h_s_1d(1:npti) * a_i_1d(1:npti) 461 sv_i_1d(1:npti) = s_i_1d(1:npti) * v_i_1d(1:npti) 462 v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti) 463 v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti) 464 sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 465 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 462 466 463 467 CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i ) 464 468 CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl) ) 465 CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,kl) )466 CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s(:,:,kl) )469 CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl) ) 470 CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl) ) 467 471 CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl) ) 468 CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i(:,:,kl) )472 CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl) ) 469 473 DO jk = 1, nlay_s 470 CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )471 CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )474 CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) ) 475 CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) ) 472 476 END DO 473 477 DO jk = 1, nlay_i 474 CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) ) 475 CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) ) 476 CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 477 END DO 478 CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl) ) 479 CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl) ) 480 CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 481 END DO 482 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) ) 483 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) ) 484 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 478 485 ! 479 486 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) … … 491 498 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 492 499 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 500 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 493 501 ! 494 502 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) … … 526 534 CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) ) 527 535 CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 536 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) ) 528 537 ! 529 538 END SELECT -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
r8565 r8623 148 148 ! Contribution to mass flux 149 149 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice 150 dh_s_mlt(ji) = dh_s_mlt(ji) - h_s_1d(ji) 150 151 ! updates 151 h_s_1d(ji) = 0._wp152 h_s_1d(ji) = 0._wp 152 153 e_s_1d (ji,1) = 0._wp 153 154 t_s_1d (ji,1) = rt0 … … 211 212 ! snow melting only = water into the ocean (then without snow precip), >0 212 213 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 214 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 213 215 ! updates available heat + precipitations after melting 214 216 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) ) … … 233 235 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 234 236 ! snow melting only = water into the ocean (then without snow precip) 235 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 237 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 238 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk) 236 239 ! updates available heat + thickness 237 zq_su (ji) 240 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 238 241 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 239 242 END DO … … 571 574 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 572 575 dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 573 h_s_1d (ji) = h_s_1d(ji) + zdeltah(ji,1)576 h_s_1d (ji) = h_s_1d(ji) + zdeltah(ji,1) 574 577 575 578 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) … … 577 580 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 578 581 ! Contribution to mass flux 579 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 582 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 583 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 580 584 ! 581 585 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90
r8598 r8623 6 6 !! history : ! Original code by Daniela Flocco and Adrian Turner 7 7 !! 1.0 ! 2012 (O. Lecomte) Adaptation for scientific tests (NEMO3.1) 8 !! 2.0 ! 2017 (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO 3.68 !! 2.0 ! 2017 (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO4 9 9 !!---------------------------------------------------------------------- 10 10 #if defined key_lim3 11 11 !!---------------------------------------------------------------------- 12 !! 'key_lim3' : LIM3sea-ice model12 !! 'key_lim3' : ESIM sea-ice model 13 13 !!---------------------------------------------------------------------- 14 14 !! ice_thd_pnd_init : some initialization and namelist read 15 15 !! ice_thd_pnd : main calling routine 16 16 !!---------------------------------------------------------------------- 17 USE phycst ! physical constants 18 USE dom_oce ! ocean space and time domain 19 USE ice ! LIM-3 variables 17 USE phycst ! physical constants 18 USE dom_oce ! ocean space and time domain 19 USE ice ! sea-ice: variables 20 USE ice1D ! sea-ice: thermodynamics variables 21 USE icetab ! sea-ice: 1D <==> 2D transformation 20 22 ! 21 USE lbclnk ! lateral boundary conditions - MPP exchanges 22 USE lib_mpp ! MPP library 23 USE in_out_manager ! I/O manager 24 USE lib_fortran ! glob_sum 25 USE timing ! Timing 23 USE in_out_manager ! I/O manager 24 USE lib_mpp ! MPP library 25 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 26 USE timing ! Timing 26 27 27 28 IMPLICIT NONE … … 40 41 # include "vectopt_loop_substitute.h90" 41 42 !!---------------------------------------------------------------------- 42 !! NEMO/ LIM3 4.0 , UCL - NEMO Consortium (2011)43 !! $Id: limdyn.F90 6994 2016-10-05 13:07:10Z clem $43 !! NEMO/ICE 4.0 , NEMO Consortium (2017) 44 !! $Id: icethd_pnd.F90 8420 2017-10-05 13:07:10Z clem $ 44 45 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 45 46 !!---------------------------------------------------------------------- 46 47 CONTAINS 47 48 48 SUBROUTINE ice_thd_pnd ( kt )49 SUBROUTINE ice_thd_pnd 49 50 !!------------------------------------------------------------------- 50 51 !! *** ROUTINE ice_thd_pnd *** … … 52 53 !! ** Purpose : change melt pond fraction 53 54 !! 54 !! ** Method : brut alforce55 !! ** Method : brut force 55 56 !! 56 57 !! ** Action : - 57 58 !! - 58 !!-------------------------------------------------------------------59 INTEGER, INTENT(in) :: kt ! number of iteration60 INTEGER :: ji, jj, jl ! dummy loop indices61 59 !!------------------------------------------------------------------- 62 60 … … 96 94 !! 97 95 !!------------------------------------------------------------------- 98 99 WHERE ( a_i(:,:,:) > 0._wp .AND. t_su(:,:,:) >= rt0 ) 100 a_ip_frac = rn_apnd 101 h_ip = rn_hpnd 102 v_ip = a_ip_frac * a_i * h_ip 103 a_ip = a_ip_frac * a_i 104 ELSEWHERE 105 a_ip = 0._wp 106 h_ip = 0._wp 107 v_ip = 0._wp 108 a_ip_frac = 0._wp 109 END WHERE 96 INTEGER :: ji ! loop indices 97 !!------------------------------------------------------------------- 98 DO ji = 1, npti 99 100 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 101 a_ip_frac_1d(ji) = rn_apnd 102 h_ip_1d(ji) = rn_hpnd 103 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 104 ELSE 105 a_ip_frac_1d(ji) = 0._wp 106 h_ip_1d(ji) = 0._wp 107 a_ip_1d(ji) = 0._wp 108 ENDIF 109 110 END DO 110 111 111 112 END SUBROUTINE pnd_CST … … 117 118 !! ** Purpose : Compute melt pond evolution 118 119 !! 119 !! ** Method : Empirical method. A fraction of meltwater is accumulated 120 !! in pond volume. It is then released exponentially when 121 !! surface is freezing. 120 !! ** Method : Empirical method. A fraction of meltwater is accumulated in ponds 121 !! and sent to ocean when surface is freezing 122 !! 123 !! pond growth: Vp = Vp + dVmelt 124 !! with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i 125 !! pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) 126 !! with Tp = -2degC 122 127 !! 123 128 !! ** Tunable parameters : (no real expertise yet, ideas?) … … 131 136 !! 132 137 !!------------------------------------------------------------------- 133 134 INTEGER, DIMENSION(jpij) :: indxi ! compressed indices for cells with ice melting135 INTEGER, DIMENSION(jpij) :: indxj !136 137 REAL(wp), DIMENSION(jpi,jpj) :: zwfx_mlw ! available meltwater for melt ponding138 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zrfrac ! fraction of available meltwater retained for melt ponding139 140 138 REAL(wp), PARAMETER :: zrmin = 0.15_wp ! minimum fraction of available meltwater retained for melt ponding 141 139 REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum '' '' '' '' '' 142 REAL(wp), PARAMETER :: zrexp = 0.01_wp ! rate constant to refreeze melt ponds143 140 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 144 145 REAL(wp) :: zhi ! dummy ice thickness 146 REAL(wp) :: z hs ! dummy snow depth147 REAL(wp) :: z Tp, z1_Tp ! reference temperature148 REAL(wp) :: z dTs ! dummy temperature difference141 REAL(wp), PARAMETER :: zTp = -2._wp ! 142 143 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 144 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding 145 REAL(wp) :: z1_Tp ! reference temperature 149 146 REAL(wp) :: z1_rhofw ! inverse freshwater density 150 147 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 151 REAL(wp) :: zvp old ! dummy pond volume152 153 INTEGER :: ji, jj, jl, ij ! loop indices 154 INTEGER :: icells ! size of dummy array148 REAL(wp) :: zvpnd ! dummy pond volume 149 REAL(wp) :: zfac, zdum 150 151 INTEGER :: ji ! loop indices 155 152 !!------------------------------------------------------------------- 156 153 z1_rhofw = 1. / rhofw 157 154 z1_zpnd_aspect = 1. / zpnd_aspect 158 zTp = -2.159 155 z1_Tp = 1._wp / zTp 160 156 161 a_ip_frac(:,:,:) = 0._wp 162 h_ip (:,:,:) = 0._wp 163 164 !------------------------------------------------------------------ 165 ! Available melt water for melt ponding and corresponding fraction 166 !------------------------------------------------------------------ 167 168 zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp ) ! available meltwater for melt ponding 169 170 ! NB: zwfx_mlw can be slightly negative for very small values (why?) 171 ! This can in some occasions give negative 172 ! v_ip in the first category, which then gives crazy pond 173 ! fractions and crashes the code as soon as the melt-pond 174 ! radiative coupling is activated 175 ! if we understand and remove why wfx_sum or wfx_snw could be 176 ! negative, then, we can remove the MAX 177 ! NB: I now changed to wfx_snw_sum, this may fix the problem. 178 ! We should check 179 180 zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:) 181 182 DO jl = 1, jpl 183 184 !------------------------------------------------------------------------------ 185 ! Identify grid cells where ponds should be updated (can probably be improved) 186 !------------------------------------------------------------------------------ 187 indxi(:) = 0 188 indxj(:) = 0 189 icells = 0 190 DO jj = 1, jpj 191 DO ji = 1, jpi 192 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 193 icells = icells + 1 194 indxi(icells) = ji 195 indxj(icells) = jj 196 ENDIF 197 END DO 198 END DO 199 200 DO ij = 1, icells 201 202 ji = indxi(ij) 203 jj = indxj(ij) 204 205 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 206 zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl) 207 208 IF ( zhi < rn_himin) THEN !--- Remove ponds on thin ice if ice is too thin 209 210 a_ip(ji,jj,jl) = 0._wp !--- Dump ponds 211 v_ip(ji,jj,jl) = 0._wp 212 a_ip_frac(ji,jj,jl) = 0._wp 213 h_ip(ji,jj,jl) = 0._wp 214 215 IF( ln_pnd_fwb ) wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) !--- Give freshwater to the ocean 216 217 ELSE !--- Update pond characteristics 218 219 !--- Add retained melt water to melt ponds 220 ! v_ip should never be positive, otherwise code crashes 221 ! MV: as far as I saw, UM5 can create very small negative v_ip values 222 ! hence I added the max, which was not required with Prather (1 yr run) 223 v_ip(ji,jj,jl) = MAX( v_ip(ji,jj,jl), 0._wp ) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice 224 225 !--- Shrink pond due to refreezing 226 zdTs = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. ) 227 228 zvpold = v_ip(ji,jj,jl) 229 230 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) * EXP( zrexp * zdTs * z1_Tp ) 231 232 !--- Dump meltwater due to refreezing ( of course this is wrong 233 !--- but this parameterization is too simple ) 234 IF ( ln_pnd_fwb ) THEN 235 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice 236 ENDIF 237 238 a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) ) 239 !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative 240 241 h_ip(ji,jj,jl) = zpnd_aspect * a_ip_frac(ji,jj,jl) 242 243 a_ip(ji,jj,jl) = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 244 157 DO ji = 1, npti 158 ! !--------------------------------! 159 IF( h_i_1d(ji) < rn_himin) THEN ! Case ice thickness < rn_himin ! 160 ! !--------------------------------! 161 !--- Remove ponds on thin ice and send water into the ocean 162 IF( ln_pnd_fwb ) THEN 163 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + h_ip_1d(ji) * a_ip_1d(ji) * rhofw * r1_rdtice 245 164 ENDIF 246 247 END DO 248 165 a_ip_1d(ji) = 0._wp 166 a_ip_frac_1d(ji) = 0._wp 167 h_ip_1d(ji) = 0._wp 168 ! !--------------------------------! 169 ELSE ! Case ice thickness >= rn_himin ! 170 ! !--------------------------------! 171 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! record pond volume at previous time step 172 zvpnd = v_ip_1d(ji) 173 174 ! available meltwater for melt ponding [m, >0] and fraction 175 zdv_mlt = -( dh_i_surf(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji) 176 zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) !!clem CICE doc 177 !!! zfr_mlt = zrmin + zrmax * a_i_1d(ji) !!clem Holland paper 178 179 !--- Pond gowth 180 ! v_ip should never be negative, otherwise code crashes 181 ! MV: as far as I saw, UM5 can create very small negative v_ip values (not Prather) 182 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 183 184 IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN 185 ! melt ponds mass flux (<0) 186 zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice 187 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 188 189 ! adjust ice/snow melting (>0) to take into account ponding 190 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 191 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 192 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 193 ENDIF 194 195 !--- Pond contraction (due to refreezing) 196 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 197 198 ! --- Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 199 ! h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i 200 a_ip_1d(ji) = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) ) 201 !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative 202 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 203 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 204 205 !! clem: there is no clever way to conserve mass here... 206 207 ! IF( ln_pnd_fwb ) THEN 208 ! ! melt ponds mass flux (>0 when ponds shrink) 209 ! zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw * r1_rdtice 210 ! wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac!! 211 ! 212 ! ! adjust ice/snow 213 ! zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw / ( a_i_1d(ji)*h_s_1d(ji)*rhosn + a_i_1d(ji)*h_i_1d(ji)*rhoic )! 214 ! 215 ! wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + h_s_1d(ji) * zfac * a_i_1d(ji) * rhosn * r1_rdtice 216 ! wfx_sum_1d(ji) = wfx_sum_1d(ji) + h_i_1d(ji) * zfac * a_i_1d(ji) * rhoic * r1_rdtice 217 ! 218 ! !!h_s_1d(ji) = h_s_1d(ji) * ( 1._wp + zfac ) 219 ! !!h_i_1d(ji) = h_i_1d(ji) * ( 1._wp + zfac )! 220 ! 221 ! ENDIF 222 223 ENDIF 249 224 END DO 250 251 !--- Remove retained meltwater from surface fluxes252 253 IF ( ln_pnd_fwb ) THEN254 255 wfx_snw_sum(:,:) = wfx_snw_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )256 wfx_sum(:,:) = wfx_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )257 258 ENDIF259 225 260 226 END SUBROUTINE pnd_H12 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
r8597 r8623 164 164 ! mass flux from ice/ocean 165 165 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 166 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 167 168 IF ( ln_pnd_fwb ) wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj) 166 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 169 167 170 168 ! add the snow melt water to snow mass flux to the ocean -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8597 r8623 255 255 !!------------------------------------------------------------------- 256 256 ! 257 v_i (:,:,:) = h_i(:,:,:) * a_i(:,:,:) 258 v_s (:,:,:) = h_s(:,:,:) * a_i(:,:,:) 259 sv_i(:,:,:) = s_i(:,:,:) * v_i(:,:,:) 257 v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:) 258 v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:) 259 sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:) 260 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 260 261 ! 261 262 END SUBROUTINE ice_var_eqv2glo … … 512 513 ! open water = 1 if at_i=0 513 514 WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp 515 514 516 ! 515 517 END SUBROUTINE ice_var_zapsmall
Note: See TracChangeset
for help on using the changeset viewer.
