Changeset 5600 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
- Timestamp:
- 2015-07-15T17:46:12+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r5034 r5600 30 30 !!====================================================================== 31 31 !! History : - ! 2006-01 (M. Vancoppenolle) Original code 32 !! 4.0! 2011-02 (G. Madec) dynamical allocation32 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 33 33 !!---------------------------------------------------------------------- 34 34 #if defined key_lim3 … … 36 36 !! 'key_lim3' LIM3 sea-ice model 37 37 !!---------------------------------------------------------------------- 38 !! lim_var_agg :39 !! lim_var_glo2eqv :40 !! lim_var_eqv2glo :41 !! lim_var_salprof :42 !! lim_var_salprof1d :43 !! lim_var_bv :44 !!----------------------------------------------------------------------45 38 USE par_oce ! ocean parameters 46 39 USE phycst ! physical constants (ocean directory) 47 40 USE sbc_oce ! Surface boundary condition: ocean fields 48 41 USE ice ! ice variables 49 USE par_ice ! ice parameters50 42 USE thd_ice ! ice variables (thermodynamics) 51 43 USE dom_ice ! ice domain … … 58 50 PRIVATE 59 51 60 PUBLIC lim_var_agg ! 61 PUBLIC lim_var_glo2eqv ! 62 PUBLIC lim_var_eqv2glo ! 63 PUBLIC lim_var_salprof ! 64 PUBLIC lim_var_icetm ! 65 PUBLIC lim_var_bv ! 66 PUBLIC lim_var_salprof1d ! 52 PUBLIC lim_var_agg 53 PUBLIC lim_var_glo2eqv 54 PUBLIC lim_var_eqv2glo 55 PUBLIC lim_var_salprof 56 PUBLIC lim_var_icetm 57 PUBLIC lim_var_bv 58 PUBLIC lim_var_salprof1d 59 PUBLIC lim_var_zapsmall 60 PUBLIC lim_var_itd 67 61 68 62 !!---------------------------------------------------------------------- 69 !! NEMO/LIM3 4.0, UCL - NEMO Consortium (2011)63 !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011) 70 64 !! $Id$ 71 65 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 129 123 DO jj = 1, jpj 130 124 DO ji = 1, jpi 131 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content132 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi 10 ) )133 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi 10 ) * rswitch ! ice salinity134 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi 10 ) )135 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi 10 ) * rswitch ! ice age125 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 126 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) ) 127 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch ! ice salinity 128 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) ) 129 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi20 ) * rswitch ! ice age 136 130 END DO 137 131 END DO … … 167 161 DO jj = 1, jpj 168 162 DO ji = 1, jpi 169 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes170 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * rswitch171 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * rswitch172 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * rswitch173 END DO 174 END DO 175 END DO 176 177 IF( n um_sal == 2 )THEN163 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes 164 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 165 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 166 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 167 END DO 168 END DO 169 END DO 170 171 IF( nn_icesal == 2 )THEN 178 172 DO jl = 1, jpl 179 173 DO jj = 1, jpj 180 174 DO ji = 1, jpi 181 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 182 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * rswitch 175 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes 176 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch 177 ! ! bounding salinity 178 sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin ) 183 179 END DO 184 180 END DO … … 191 187 ! Ice temperatures 192 188 !------------------- 193 !CDIR NOVERRCHK 194 DO jl = 1, jpl 195 !CDIR NOVERRCHK 189 DO jl = 1, jpl 196 190 DO jk = 1, nlay_i 197 !CDIR NOVERRCHK198 191 DO jj = 1, jpj 199 !CDIR NOVERRCHK200 192 DO ji = 1, jpi 201 193 ! ! Energy of melting q(S,T) [J.m-3] 202 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) ! rswitch = 0 if no ice and 1 if yes 203 zq_i = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 204 zq_i = zq_i * unit_fac !convert units 205 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt ! Ice layer melt temperature 194 rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 195 zq_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 196 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0 ! Ice layer melt temperature 206 197 ! 207 198 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 208 zbbb = ( rcp - cpic ) * ( ztmelts - rt t ) + zq_i /rhoic - lfus209 zccc = lfus * (ztmelts-rt t)199 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus 200 zccc = lfus * (ztmelts-rt0) 210 201 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 211 t_i(ji,jj,jk,jl) = rt t+ rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )212 t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt202 t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 203 t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) ) ! -100 < t_i < ztmelts 213 204 END DO 214 205 END DO … … 226 217 DO ji = 1, jpi 227 218 !Energy of melting q(S,T) [J.m-3] 228 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) ) ! rswitch = 0 if no ice and 1 if yes 229 zq_s = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 230 zq_s = zq_s * unit_fac ! convert units 219 rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 220 zq_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 231 221 ! 232 t_s(ji,jj,jk,jl) = rt t+ rswitch * ( - zfac1 * zq_s + zfac2 )233 t_s(ji,jj,jk,jl) = MIN( rt t, MAX( 173.15, t_s(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt222 t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 ) 223 t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) ) ! -100 < t_s < rt0 234 224 END DO 235 225 END DO … … 240 230 ! Mean temperature 241 231 !------------------- 232 vt_i (:,:) = 0._wp 233 DO jl = 1, jpl 234 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 235 END DO 236 242 237 tm_i(:,:) = 0._wp 243 238 DO jl = 1, jpl … … 245 240 DO jj = 1, jpj 246 241 DO ji = 1, jpi 247 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 248 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 249 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 250 END DO 251 END DO 252 END DO 253 END DO 242 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 243 tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) & 244 & / MAX( vt_i(ji,jj) , epsi10 ) 245 END DO 246 END DO 247 END DO 248 END DO 249 tm_i = tm_i + rt0 254 250 ! 255 251 END SUBROUTINE lim_var_glo2eqv … … 270 266 v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:) 271 267 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 272 oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:)273 268 ! 274 269 END SUBROUTINE lim_var_eqv2glo … … 281 276 !! ** Purpose : computes salinity profile in function of bulk salinity 282 277 !! 283 !! ** Method : If bulk salinity greater than s_i_1,278 !! ** Method : If bulk salinity greater than zsi1, 284 279 !! the profile is assumed to be constant (S_inf) 285 !! If bulk salinity lower than s_i_0,280 !! If bulk salinity lower than zsi0, 286 281 !! the profile is linear with 0 at the surface (S_zero) 287 !! If it is between s_i_0 and s_i_1, it is a282 !! If it is between zsi0 and zsi1, it is a 288 283 !! alpha-weighted linear combination of s_inf and s_zero 289 284 !! 290 !! ** References : Vancoppenolle et al., 2007 (in preparation)285 !! ** References : Vancoppenolle et al., 2007 291 286 !!------------------------------------------------------------------ 292 287 INTEGER :: ji, jj, jk, jl ! dummy loop index 293 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac, zsal ! local scalar 294 REAL(wp) :: zswi0, zswi01, zswibal, zargtemp , zs_zero ! - - 295 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha ! 3D pointer 288 REAL(wp) :: zfac0, zfac1, zsal 289 REAL(wp) :: zswi0, zswi01, zargtemp , zs_zero 290 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha 291 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 292 REAL(wp), PARAMETER :: zsi1 = 4.5_wp 296 293 !!------------------------------------------------------------------ 297 294 … … 301 298 ! Vertically constant, constant in time 302 299 !--------------------------------------- 303 IF( n um_sal == 1 ) s_i(:,:,:,:) = bulk_sal300 IF( nn_icesal == 1 ) s_i(:,:,:,:) = rn_icesal 304 301 305 302 !----------------------------------- 306 303 ! Salinity profile, varying in time 307 304 !----------------------------------- 308 IF( n um_sal == 2 ) THEN305 IF( nn_icesal == 2 ) THEN 309 306 ! 310 307 DO jk = 1, nlay_i … … 315 312 DO jj = 1, jpj 316 313 DO ji = 1, jpi 317 z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) ) 318 END DO 319 END DO 320 END DO 321 ! 322 dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 ) ! Weighting factor between zs_zero and zs_inf 323 dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 ) 314 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) ) 315 z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) ) 316 END DO 317 END DO 318 END DO 319 ! 320 zfac0 = 1._wp / ( zsi0 - zsi1 ) ! Weighting factor between zs_zero and zs_inf 321 zfac1 = zsi1 / ( zsi1 - zsi0 ) 324 322 ! 325 323 zalpha(:,:,:) = 0._wp … … 327 325 DO jj = 1, jpj 328 326 DO ji = 1, jpi 329 ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise330 zswi0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i(ji,jj,jl) ) )331 ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws332 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i(ji,jj,jl) ) )333 ! If 2.sm_i GE sss_m then zswibal= 1327 ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise 328 zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i(ji,jj,jl) ) ) 329 ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 330 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i(ji,jj,jl) ) ) 331 ! If 2.sm_i GE sss_m then rswitch = 1 334 332 ! this is to force a constant salinity profile in the Baltic Sea 335 zswibal= MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )336 zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )337 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal)338 END DO 339 END DO 340 END DO 341 342 dummy_fac = 1._wp / REAL( nlay_i )! Computation of the profile333 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 334 zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 ) 335 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch ) 336 END DO 337 END DO 338 END DO 339 340 ! Computation of the profile 343 341 DO jl = 1, jpl 344 342 DO jk = 1, nlay_i … … 346 344 DO ji = 1, jpi 347 345 ! ! linear profile with 0 at the surface 348 zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac346 zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i 349 347 ! ! weighting the profile 350 348 s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 351 END DO ! ji 352 END DO ! jj 353 END DO ! jk 354 END DO ! jl 355 ! 356 ENDIF ! num_sal 349 ! ! bounding salinity 350 s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) ) 351 END DO 352 END DO 353 END DO 354 END DO 355 ! 356 ENDIF ! nn_icesal 357 357 358 358 !------------------------------------------------------- … … 360 360 !------------------------------------------------------- 361 361 362 IF( n um_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)362 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 363 363 ! 364 364 sm_i(:,:,:) = 2.30_wp 365 365 ! 366 366 DO jl = 1, jpl 367 !CDIR NOVERRCHK368 367 DO jk = 1, nlay_i 369 zargtemp = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)368 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 370 369 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 371 370 s_i(:,:,jk,jl) = zsal … … 373 372 END DO 374 373 ! 375 ENDIF ! n um_sal374 ENDIF ! nn_icesal 376 375 ! 377 376 CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha ) … … 390 389 391 390 ! Mean sea ice temperature 391 vt_i (:,:) = 0._wp 392 DO jl = 1, jpl 393 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 394 END DO 395 392 396 tm_i(:,:) = 0._wp 393 397 DO jl = 1, jpl … … 395 399 DO jj = 1, jpj 396 400 DO ji = 1, jpi 397 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 398 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 399 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 400 END DO 401 END DO 402 END DO 403 END DO 401 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 402 tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) & 403 & / MAX( vt_i(ji,jj) , epsi10 ) 404 END DO 405 END DO 406 END DO 407 END DO 408 tm_i = tm_i + rt0 404 409 405 410 END SUBROUTINE lim_var_icetm … … 420 425 !!------------------------------------------------------------------ 421 426 ! 427 vt_i (:,:) = 0._wp 428 DO jl = 1, jpl 429 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 430 END DO 431 422 432 bv_i(:,:) = 0._wp 423 433 DO jl = 1, jpl … … 425 435 DO jj = 1, jpj 426 436 DO ji = 1, jpi 427 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt t) + epsi10 ) ) )428 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt t, - epsi10 ) &429 & * v_i(ji,jj,jl) / REAL(nlay_i,wp)430 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi 10 ) ) )431 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi 10 )437 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) ) 438 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) & 439 & * v_i(ji,jj,jl) * r1_nlay_i 440 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) ) ) 441 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi20 ) 432 442 END DO 433 443 END DO … … 448 458 ! 449 459 INTEGER :: ji, jk ! dummy loop indices 450 INTEGER :: ii, ij ! local integers451 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal ! local scalars452 REAL(wp) :: zalpha, zswi0, zswi01, zs wibal, zs_zero ! - -460 INTEGER :: ii, ij ! local integers 461 REAL(wp) :: zfac0, zfac1, zargtemp, zsal ! local scalars 462 REAL(wp) :: zalpha, zswi0, zswi01, zs_zero ! - - 453 463 ! 454 464 REAL(wp), POINTER, DIMENSION(:) :: z_slope_s 465 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 466 REAL(wp), PARAMETER :: zsi1 = 4.5_wp 455 467 !!--------------------------------------------------------------------- 456 468 … … 460 472 ! Vertically constant, constant in time 461 473 !--------------------------------------- 462 IF( n um_sal == 1 ) s_i_1d(:,:) = bulk_sal474 IF( nn_icesal == 1 ) s_i_1d(:,:) = rn_icesal 463 475 464 476 !------------------------------------------------------ … … 466 478 !------------------------------------------------------ 467 479 468 IF( n um_sal == 2 ) THEN480 IF( nn_icesal == 2 ) THEN 469 481 ! 470 482 DO ji = kideb, kiut ! Slope of the linear profile zs_zero 471 z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) ) 483 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 484 z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) ) 472 485 END DO 473 486 474 487 ! Weighting factor between zs_zero and zs_inf 475 488 !--------------------------------------------- 476 dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 ) 477 dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 ) 478 dummy_fac2 = 1._wp / REAL(nlay_i,wp) 479 480 !CDIR NOVERRCHK 489 zfac0 = 1._wp / ( zsi0 - zsi1 ) 490 zfac1 = zsi1 / ( zsi1 - zsi0 ) 481 491 DO jk = 1, nlay_i 482 !CDIR NOVERRCHK483 492 DO ji = kideb, kiut 484 493 ii = MOD( npb(ji) - 1 , jpi ) + 1 485 494 ij = ( npb(ji) - 1 ) / jpi + 1 486 ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise487 zswi0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_1d(ji) ) )488 ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws489 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_1d(ji) ) )490 ! if 2.sm_i GE sss_m then zswibal= 1495 ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise 496 zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i_1d(ji) ) ) 497 ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 498 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 499 ! if 2.sm_i GE sss_m then rswitch = 1 491 500 ! this is to force a constant salinity profile in the Baltic Sea 492 zswibal= MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )501 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 493 502 ! 494 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zswibal)503 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 ) ) * ( 1._wp - rswitch ) 495 504 ! 496 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2505 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 497 506 ! weighting the profile 498 507 s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 499 END DO ! ji 500 END DO ! jk 501 502 ENDIF ! num_sal 508 ! bounding salinity 509 s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) ) 510 END DO 511 END DO 512 513 ENDIF 503 514 504 515 !------------------------------------------------------- … … 506 517 !------------------------------------------------------- 507 518 508 IF( n um_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)519 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 509 520 ! 510 521 sm_i_1d(:) = 2.30_wp 511 522 ! 512 !CDIR NOVERRCHK513 523 DO jk = 1, nlay_i 514 zargtemp = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)515 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ))524 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 525 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) ) 516 526 DO ji = kideb, kiut 517 527 s_i_1d(ji,jk) = zsal … … 524 534 ! 525 535 END SUBROUTINE lim_var_salprof1d 536 537 SUBROUTINE lim_var_zapsmall 538 !!------------------------------------------------------------------- 539 !! *** ROUTINE lim_var_zapsmall *** 540 !! 541 !! ** Purpose : Remove too small sea ice areas and correct fluxes 542 !! 543 !! history : LIM3.5 - 01-2014 (C. Rousset) original code 544 !!------------------------------------------------------------------- 545 INTEGER :: ji, jj, jl, jk ! dummy loop indices 546 REAL(wp) :: zsal, zvi, zvs, zei, zes 547 !!------------------------------------------------------------------- 548 at_i (:,:) = 0._wp 549 DO jl = 1, jpl 550 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 551 END DO 552 553 DO jl = 1, jpl 554 555 !----------------------------------------------------------------- 556 ! Zap ice energy and use ocean heat to melt ice 557 !----------------------------------------------------------------- 558 DO jk = 1, nlay_i 559 DO jj = 1 , jpj 560 DO ji = 1 , jpi 561 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 562 rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj ) - epsi10 ) ) * rswitch 563 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 564 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch & 565 & / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 566 zei = e_i(ji,jj,jk,jl) 567 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch 568 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 569 ! update exchanges with ocean 570 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 571 END DO 572 END DO 573 END DO 574 575 DO jj = 1 , jpj 576 DO ji = 1 , jpi 577 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 578 rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj ) - epsi10 ) ) * rswitch 579 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 580 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch & 581 & / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 582 zsal = smv_i(ji,jj, jl) 583 zvi = v_i (ji,jj, jl) 584 zvs = v_s (ji,jj, jl) 585 zes = e_s (ji,jj,1,jl) 586 !----------------------------------------------------------------- 587 ! Zap snow energy 588 !----------------------------------------------------------------- 589 t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 590 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch 591 592 !----------------------------------------------------------------- 593 ! zap ice and snow volume, add water and salt to ocean 594 !----------------------------------------------------------------- 595 ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) 596 a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch 597 v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch 598 v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch 599 t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) 600 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch 601 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 602 603 ! update exchanges with ocean 604 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 605 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 606 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 607 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 608 END DO 609 END DO 610 END DO 611 612 ! to be sure that at_i is the sum of a_i(jl) 613 at_i (:,:) = 0._wp 614 DO jl = 1, jpl 615 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 616 END DO 617 618 ! open water = 1 if at_i=0 619 DO jj = 1, jpj 620 DO ji = 1, jpi 621 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 622 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 623 END DO 624 END DO 625 626 ! 627 END SUBROUTINE lim_var_zapsmall 628 629 SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i ) 630 !!------------------------------------------------------------------ 631 !! *** ROUTINE lim_var_itd *** 632 !! 633 !! ** Purpose : converting 1-cat ice to multiple ice categories 634 !! 635 !! ice thickness distribution follows a gaussian law 636 !! around the concentration of the most likely ice thickness 637 !! (similar as limistate.F90) 638 !! 639 !! ** Method: Iterative procedure 640 !! 641 !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 642 !! 643 !! 2) Check whether the distribution conserves area and volume, positivity and 644 !! category boundaries 645 !! 646 !! 3) If not (input ice is too thin), the last category is empty and 647 !! the number of categories is reduced (jpl-1) 648 !! 649 !! 4) Iterate until ok (SUM(itest(:) = 4) 650 !! 651 !! ** Arguments : zhti: 1-cat ice thickness 652 !! zhts: 1-cat snow depth 653 !! zai : 1-cat ice concentration 654 !! 655 !! ** Output : jpl-cat 656 !! 657 !! (Example of application: BDY forcings when input are cell averaged) 658 !! 659 !!------------------------------------------------------------------- 660 !! History : LIM3.5 - 2012 (M. Vancoppenolle) Original code 661 !! 2014 (C. Rousset) Rewriting 662 !!------------------------------------------------------------------- 663 !! Local variables 664 INTEGER :: ji, jk, jl ! dummy loop indices 665 INTEGER :: ijpij, i_fill, jl0 666 REAL(wp) :: zarg, zV, zconv, zdh 667 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables 668 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables 669 INTEGER , POINTER, DIMENSION(:) :: itest 670 671 CALL wrk_alloc( 4, itest ) 672 !-------------------------------------------------------------------- 673 ! initialisation of variables 674 !-------------------------------------------------------------------- 675 ijpij = SIZE(zhti,1) 676 zht_i(1:ijpij,1:jpl) = 0._wp 677 zht_s(1:ijpij,1:jpl) = 0._wp 678 za_i (1:ijpij,1:jpl) = 0._wp 679 680 ! ---------------------------------------- 681 ! distribution over the jpl ice categories 682 ! ---------------------------------------- 683 DO ji = 1, ijpij 684 685 IF( zhti(ji) > 0._wp ) THEN 686 687 ! initialisation of tests 688 itest(:) = 0 689 690 i_fill = jpl + 1 !==================================== 691 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 692 ! iteration !==================================== 693 i_fill = i_fill - 1 694 695 ! initialisation of ice variables for each try 696 zht_i(ji,1:jpl) = 0._wp 697 za_i (ji,1:jpl) = 0._wp 698 699 ! *** case very thin ice: fill only category 1 700 IF ( i_fill == 1 ) THEN 701 zht_i(ji,1) = zhti(ji) 702 za_i (ji,1) = zai (ji) 703 704 ! *** case ice is thicker: fill categories >1 705 ELSE 706 707 ! Fill ice thicknesses except the last one (i_fill) by hmean 708 DO jl = 1, i_fill - 1 709 zht_i(ji,jl) = hi_mean(jl) 710 END DO 711 712 ! find which category (jl0) the input ice thickness falls into 713 jl0 = i_fill 714 DO jl = 1, i_fill 715 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 716 jl0 = jl 717 CYCLE 718 ENDIF 719 END DO 720 721 ! Concentrations in the (i_fill-1) categories 722 za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 723 DO jl = 1, i_fill - 1 724 IF ( jl == jl0 ) CYCLE 725 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 726 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 727 END DO 728 729 ! Concentration in the last (i_fill) category 730 za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 731 732 ! Ice thickness in the last (i_fill) category 733 zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 734 zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 735 736 ENDIF ! case ice is thick or thin 737 738 !--------------------- 739 ! Compatibility tests 740 !--------------------- 741 ! Test 1: area conservation 742 zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 743 IF ( zconv < epsi06 ) itest(1) = 1 744 745 ! Test 2: volume conservation 746 zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 747 IF ( zconv < epsi06 ) itest(2) = 1 748 749 ! Test 3: thickness of the last category is in-bounds ? 750 IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 751 752 ! Test 4: positivity of ice concentrations 753 itest(4) = 1 754 DO jl = 1, i_fill 755 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 756 END DO 757 !============================ 758 END DO ! end iteration on categories 759 !============================ 760 ENDIF ! if zhti > 0 761 END DO ! i loop 762 763 ! ------------------------------------------------ 764 ! Adding Snow in each category where za_i is not 0 765 ! ------------------------------------------------ 766 DO jl = 1, jpl 767 DO ji = 1, ijpij 768 IF( za_i(ji,jl) > 0._wp ) THEN 769 zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 770 ! In case snow load is in excess that would lead to transformation from snow to ice 771 ! Then, transfer the snow excess into the ice (different from limthd_dh) 772 zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 773 ! recompute ht_i, ht_s avoiding out of bounds values 774 zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) 775 zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn ) 776 ENDIF 777 ENDDO 778 ENDDO 779 780 CALL wrk_dealloc( 4, itest ) 781 ! 782 END SUBROUTINE lim_var_itd 783 526 784 527 785 #else … … 542 800 SUBROUTINE lim_var_salprof1d ! Emtpy routines 543 801 END SUBROUTINE lim_var_salprof1d 802 SUBROUTINE lim_var_zapsmall 803 END SUBROUTINE lim_var_zapsmall 804 SUBROUTINE lim_var_itd 805 END SUBROUTINE lim_var_itd 544 806 #endif 545 807
Note: See TracChangeset
for help on using the changeset viewer.
