Changeset 10415
- Timestamp:
- 2018-12-19T12:28:25+01:00 (4 years ago)
- Location:
- NEMO/trunk/src
- Files:
-
- 7 edited
-
ICE/ice.F90 (modified) (2 diffs)
-
ICE/icectl.F90 (modified) (1 diff)
-
ICE/icedyn.F90 (modified) (9 diffs)
-
ICE/icedyn_rhg_evp.F90 (modified) (3 diffs)
-
ICE/icestp.F90 (modified) (2 diffs)
-
ICE/icevar.F90 (modified) (2 diffs)
-
OCE/stpctl.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/ice.F90
r10413 r10415 332 332 !! * Old values of global variables 333 333 !!---------------------------------------------------------------------- 334 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_b, sv_i_b, oa_i_b !:336 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s_b !: snow heat content337 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i_b !: ice temperatures338 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice_b, v_ice_b !: ice velocity339 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i_b !: ice concentration (total)334 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/thickness 335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_b, sv_i_b, oa_i_b !: 336 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s_b !: snow heat content 337 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i_b !: ice temperatures 338 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice_b, v_ice_b !: ice velocity 339 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i_b !: ice concentration (total) 340 340 341 341 !!---------------------------------------------------------------------- … … 441 441 ! * Old values of global variables 442 442 ii = ii + 1 443 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) ,&444 & 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) , &443 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), & 444 & 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) , & 445 445 & oa_i_b(jpi,jpj,jpl) , STAT=ierr(ii) ) 446 446 -
NEMO/trunk/src/ICE/icectl.F90
r10069 r10415 190 190 191 191 ! heat flux 192 zhfx = glob_sum( ( qt_atm_oi - qt_oce_ai - diag_heat - diag_trp_ei - diag_trp_es & 193 ! & - SUM( qevap_ice * a_i_b, dim=3 ) & !!clem: I think this line must be commented (but need check) 194 & ) * e1e2t ) * zconv 192 zhfx = glob_sum( ( qt_atm_oi - qt_oce_ai - diag_heat ) * e1e2t ) * zconv 195 193 196 194 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) -
NEMO/trunk/src/ICE/icedyn.F90
r10413 r10415 76 76 INTEGER :: ji, jj, jl ! dummy loop indices 77 77 REAL(wp) :: zcoefu, zcoefv 78 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max 78 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 79 79 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 80 80 !!-------------------------------------------------------------------- … … 101 101 DO jj = 2, jpjm1 102 102 DO ji = fs_2, fs_jpim1 103 zhi_max(ji,jj,jl) = MAX( epsi20, h_i_b(ji,jj,jl), h_i_b(ji+1,jj ,jl), h_i_b(ji ,jj+1,jl), & 104 & h_i_b(ji-1,jj ,jl), h_i_b(ji ,jj-1,jl), & 105 & h_i_b(ji+1,jj+1,jl), h_i_b(ji-1,jj-1,jl), & 106 & h_i_b(ji+1,jj-1,jl), h_i_b(ji-1,jj+1,jl) ) 107 zhs_max(ji,jj,jl) = MAX( epsi20, h_s_b(ji,jj,jl), h_s_b(ji+1,jj ,jl), h_s_b(ji ,jj+1,jl), & 108 & h_s_b(ji-1,jj ,jl), h_s_b(ji ,jj-1,jl), & 109 & h_s_b(ji+1,jj+1,jl), h_s_b(ji-1,jj-1,jl), & 110 & h_s_b(ji+1,jj-1,jl), h_s_b(ji-1,jj+1,jl) ) 103 zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj ,jl), h_ip_b(ji ,jj+1,jl), & 104 & h_ip_b(ji-1,jj ,jl), h_ip_b(ji ,jj-1,jl), & 105 & h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 106 & h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 107 zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj ,jl), h_i_b (ji ,jj+1,jl), & 108 & h_i_b (ji-1,jj ,jl), h_i_b (ji ,jj-1,jl), & 109 & h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 110 & h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 111 zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj ,jl), h_s_b (ji ,jj+1,jl), & 112 & h_s_b (ji-1,jj ,jl), h_s_b (ji ,jj-1,jl), & 113 & h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 114 & h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 111 115 END DO 112 116 END DO 113 117 END DO 114 CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1. )118 CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 115 119 ! 116 120 ! … … 118 122 119 123 CASE ( np_dynALL ) !== all dynamical processes ==! 120 CALL ice_dyn_rhg ( kt ) ! -- rheology121 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max ) ! -- advection of ice + correction on ice thickness122 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting123 CALL ice_cor ( kt , 1 ) ! -- Corrections124 CALL ice_dyn_rhg ( kt ) ! -- rheology 125 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 126 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting 127 CALL ice_cor ( kt , 1 ) ! -- Corrections 124 128 125 129 CASE ( np_dynRHGADV ) !== no ridge/raft & no corrections ==! 126 CALL ice_dyn_rhg ( kt ) ! -- rheology127 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max ) ! -- advection of ice + correction on ice thickness128 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting)130 CALL ice_dyn_rhg ( kt ) ! -- rheology 131 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 132 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 129 133 130 134 CASE ( np_dynADV1D ) !== pure advection ==! (1D) … … 163 167 !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 164 168 ! --- 165 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max ) ! -- advection of ice + correction on ice thickness169 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 166 170 167 171 ! diagnostics: divergence at T points … … 186 190 187 191 188 SUBROUTINE Hbig( phi_max, phs_max )192 SUBROUTINE Hbig( phi_max, phs_max, phip_max ) 189 193 !!------------------------------------------------------------------- 190 194 !! *** ROUTINE Hbig *** … … 203 207 !! ** input : Max thickness of the surrounding 9-points 204 208 !!------------------------------------------------------------------- 205 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi_max, phs_max ! max ice thick from surrounding 9-pts209 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 206 210 ! 207 211 INTEGER :: ji, jj, jl ! dummy loop indices 208 REAL(wp) :: zhi , zhs, zvs_excess, zfra212 REAL(wp) :: zhip, zhi, zhs, zvs_excess, zfra 209 213 !!------------------------------------------------------------------- 210 214 ! … … 216 220 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 217 221 ! 222 ! ! -- check h_ip -- ! 223 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 224 IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN 225 zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) ) 226 IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN 227 a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl) 228 ENDIF 229 ENDIF 230 ! 218 231 ! ! -- check h_i -- ! 219 232 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 220 233 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 221 !!clem zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)222 !!clem IF ( ( zdv > 0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &223 !!clem & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN224 234 IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 225 235 a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) … … 233 243 ! 234 244 wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 235 hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0245 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 236 246 ! 237 e_s(ji,jj,1 ,jl) = e_s(ji,jj,1,jl) * zfra238 v_s(ji,jj,jl) = a_i(ji,jj,jl) * phs_max(ji,jj,jl)247 e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 248 v_s(ji,jj,jl) = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 239 249 ENDIF 240 250 ! … … 242 252 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 243 253 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 244 ! this imposed mini can artificially make the snow thickness very high(if concentration decreases drastically)254 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 245 255 zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 246 256 IF( zvs_excess > 0._wp ) THEN 247 257 zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 ) 248 258 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 249 hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0259 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 250 260 ! 251 e_s(ji,jj,1 ,jl) = e_s(ji,jj,1,jl) * zfra252 v_s(ji,jj,jl) = v_s(ji,jj,jl) - zvs_excess261 e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 262 v_s(ji,jj,jl) = v_s(ji,jj,jl) - zvs_excess 253 263 ENDIF 254 264 -
NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
r10413 r10415 147 147 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 148 148 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! check convergence 149 REAL(wp), DIMENSION(jpi,jpj) :: zssh _lead_m! array used for the calculation of ice surface slope:149 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 150 150 ! ! ocean surface (ssh_m) if ice is not embedded 151 ! ! ice topsurface if ice is embedded151 ! ! ice bottom surface if ice is embedded 152 152 REAL(wp), DIMENSION(jpi,jpj) :: zCorx, zCory ! Coriolis stress array 153 153 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! Ocean-to-ice stress array … … 268 268 ! 2) Wind / ocean stress, mass terms, coriolis terms 269 269 !------------------------------------------------------------------------------! 270 271 ! == embedded sea ice: compute representative ice top surface ==!272 ! == non-embedded sea ice: use ocean surface for slope calculation ==!273 zssh _lead_m(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)270 ! sea surface height 271 ! embedded sea ice: compute representative ice top surface 272 ! non-embedded sea ice: use ocean surface for slope calculation 273 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 274 274 275 275 DO jj = 2, jpjm1 … … 309 309 310 310 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 311 zspgU(ji,jj) = - zmassU * grav * ( zssh _lead_m(ji+1,jj) - zssh_lead_m(ji,jj) ) * r1_e1u(ji,jj)312 zspgV(ji,jj) = - zmassV * grav * ( zssh _lead_m(ji,jj+1) - zssh_lead_m(ji,jj) ) * r1_e2v(ji,jj)311 zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 312 zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 313 313 314 314 ! masks -
NEMO/trunk/src/ICE/icestp.F90
r10069 r10415 194 194 CALL ice_var_agg( 2 ) ! necessary calls (at least for coupling) 195 195 ! 196 !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work197 !! moreover it should only be called at the update frequency (as in agrif_ice_update.F90)198 !# if defined key_agrif199 ! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid()200 !# endif201 196 CALL ice_update_flx( kt ) ! -- Update ocean surface mass, heat and salt fluxes 202 !# if defined key_agrif 203 ! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid() 204 !# endif 197 ! 205 198 IF( ln_icediahsb ) CALL ice_dia( kt ) ! -- Diagnostics outputs 206 199 ! … … 368 361 e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy 369 362 WHERE( a_i_b(:,:,:) >= epsi20 ) 370 h_i_b(:,:,:) = v_i_b 371 h_s_b(:,:,:) = v_s_b 363 h_i_b(:,:,:) = v_i_b(:,:,:) / a_i_b(:,:,:) ! ice thickness 364 h_s_b(:,:,:) = v_s_b(:,:,:) / a_i_b(:,:,:) ! snw thickness 372 365 ELSEWHERE 373 366 h_i_b(:,:,:) = 0._wp 374 367 h_s_b(:,:,:) = 0._wp 368 END WHERE 369 370 WHERE( a_ip(:,:,:) >= epsi20 ) 371 h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) ! ice pond thickness 372 ELSEWHERE 373 h_ip_b(:,:,:) = 0._wp 375 374 END WHERE 376 375 ! -
NEMO/trunk/src/ICE/icevar.F90
r10413 r10415 954 954 FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 955 955 !!--------------------------------------------------------------------- 956 !! *** ROUTINE rhg_evp_rst***956 !! *** ROUTINE ice_var_sshdyn *** 957 957 !! 958 958 !! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded … … 963 963 !! Sea ice-ocean coupling using a rescaled vertical coordinate z*, 964 964 !! Ocean Modelling, Volume 24, Issues 1-2, 2008 965 !!966 965 !!---------------------------------------------------------------------- 967 966 ! -
NEMO/trunk/src/OCE/stpctl.F90
r10364 r10415 50 50 !! - Print it each 50 time steps 51 51 !! - Stop the run IF problem encountered by setting indic=-3 52 !! Problems checked: |ssh| maximum larger than 10 m52 !! Problems checked: |ssh| maximum larger than 20 m 53 53 !! |U| maximum larger than 10 m/s 54 54 !! negative sea surface salinity … … 136 136 ENDIF 137 137 ! 138 IF ( zmax(1) > 15._wp .OR. & ! too large sea surface height ( > 15m )138 IF ( zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 139 139 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 140 140 & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity … … 159 159 IF(lwp) THEN 160 160 WRITE(numout,cform_err) 161 WRITE(numout,*) ' stp_ctl: |ssh| > 10 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests'161 WRITE(numout,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 162 162 WRITE(numout,*) ' ======= ' 163 163 WRITE(numout,9100) kt, zmax(1), iih , ijh
Note: See TracChangeset
for help on using the changeset viewer.
