- Timestamp:
- 2015-07-10T13:28:53+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r4688 r5581 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 !67 68 REAL(wp) :: epsi10 = 1.e-10_wp ! - -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 69 61 70 62 !!---------------------------------------------------------------------- 71 !! NEMO/LIM3 4.0, UCL - NEMO Consortium (2011)63 !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011) 72 64 !! $Id$ 73 65 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 92 84 ! 93 85 INTEGER :: ji, jj, jk, jl ! dummy loop indices 94 REAL(wp) :: zinda, zindb95 86 !!------------------------------------------------------------------ 96 87 … … 111 102 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 112 103 ! 113 zinda= MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )114 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda! ice thickness104 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 105 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch ! ice thickness 115 106 END DO 116 107 END DO … … 132 123 DO jj = 1, jpj 133 124 DO ji = 1, jpi 134 zinda = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )135 zindb = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )136 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content137 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda ! ice salinity138 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi 10 ) * zindb! 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 139 130 END DO 140 131 END DO … … 161 152 INTEGER :: ji, jj, jk, jl ! dummy loop indices 162 153 REAL(wp) :: zq_i, zaaa, zbbb, zccc, zdiscrim ! local scalars 163 REAL(wp) :: ztmelts, z indb, zq_s, zfac1, zfac2 ! - -154 REAL(wp) :: ztmelts, zq_s, zfac1, zfac2 ! - - 164 155 !!------------------------------------------------------------------ 165 156 … … 170 161 DO jj = 1, jpj 171 162 DO ji = 1, jpi 172 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes173 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * zindb174 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * zindb175 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi 10 ) * zindb176 END DO 177 END DO 178 END DO 179 180 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 181 172 DO jl = 1, jpl 182 173 DO jj = 1, jpj 183 174 DO ji = 1, jpi 184 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 185 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * zindb 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 ) 186 179 END DO 187 180 END DO … … 194 187 ! Ice temperatures 195 188 !------------------- 196 !CDIR NOVERRCHK 197 DO jl = 1, jpl 198 !CDIR NOVERRCHK 189 DO jl = 1, jpl 199 190 DO jk = 1, nlay_i 200 !CDIR NOVERRCHK201 191 DO jj = 1, jpj 202 !CDIR NOVERRCHK203 192 DO ji = 1, jpi 204 193 ! ! Energy of melting q(S,T) [J.m-3] 205 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) ! zindb = 0 if no ice and 1 if yes 206 zq_i = zindb * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 207 zq_i = zq_i * unit_fac !convert units 208 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 209 197 ! 210 198 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 211 zbbb = ( rcp - cpic ) * ( ztmelts - rt t ) + zq_i /rhoic - lfus212 zccc = lfus * (ztmelts-rt t)199 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus 200 zccc = lfus * (ztmelts-rt0) 213 201 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 214 t_i(ji,jj,jk,jl) = rt t + zindb*( - zbbb - zdiscrim ) / ( 2.0 *zaaa )215 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 216 204 END DO 217 205 END DO … … 229 217 DO ji = 1, jpi 230 218 !Energy of melting q(S,T) [J.m-3] 231 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) ) ! zindb = 0 if no ice and 1 if yes 232 zq_s = zindb * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 233 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) 234 221 ! 235 t_s(ji,jj,jk,jl) = rt t + zindb* ( - zfac1 * zq_s + zfac2 )236 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 237 224 END DO 238 225 END DO … … 243 230 ! Mean temperature 244 231 !------------------- 232 vt_i (:,:) = 0._wp 233 DO jl = 1, jpl 234 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 235 END DO 236 245 237 tm_i(:,:) = 0._wp 246 238 DO jl = 1, jpl … … 248 240 DO jj = 1, jpj 249 241 DO ji = 1, jpi 250 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 251 tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 252 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 253 END DO 254 END DO 255 END DO 256 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 257 250 ! 258 251 END SUBROUTINE lim_var_glo2eqv … … 273 266 v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:) 274 267 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 275 oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:)276 268 ! 277 269 END SUBROUTINE lim_var_eqv2glo … … 284 276 !! ** Purpose : computes salinity profile in function of bulk salinity 285 277 !! 286 !! ** Method : If bulk salinity greater than s_i_1,278 !! ** Method : If bulk salinity greater than zsi1, 287 279 !! the profile is assumed to be constant (S_inf) 288 !! If bulk salinity lower than s_i_0,280 !! If bulk salinity lower than zsi0, 289 281 !! the profile is linear with 0 at the surface (S_zero) 290 !! If it is between s_i_0 and s_i_1, it is a282 !! If it is between zsi0 and zsi1, it is a 291 283 !! alpha-weighted linear combination of s_inf and s_zero 292 284 !! 293 !! ** References : Vancoppenolle et al., 2007 (in preparation)285 !! ** References : Vancoppenolle et al., 2007 294 286 !!------------------------------------------------------------------ 295 287 INTEGER :: ji, jj, jk, jl ! dummy loop index 296 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac, zsal ! local scalar 297 REAL(wp) :: zind0, zind01, zindbal, zargtemp , zs_zero ! - - 298 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 299 293 !!------------------------------------------------------------------ 300 294 … … 304 298 ! Vertically constant, constant in time 305 299 !--------------------------------------- 306 IF( n um_sal == 1 ) s_i(:,:,:,:) = bulk_sal300 IF( nn_icesal == 1 ) s_i(:,:,:,:) = rn_icesal 307 301 308 302 !----------------------------------- 309 303 ! Salinity profile, varying in time 310 304 !----------------------------------- 311 IF( n um_sal == 2 ) THEN305 IF( nn_icesal == 2 ) THEN 312 306 ! 313 307 DO jk = 1, nlay_i … … 318 312 DO jj = 1, jpj 319 313 DO ji = 1, jpi 320 z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) ) 321 END DO 322 END DO 323 END DO 324 ! 325 dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 ) ! Weighting factor between zs_zero and zs_inf 326 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 ) 327 322 ! 328 323 zalpha(:,:,:) = 0._wp … … 330 325 DO jj = 1, jpj 331 326 DO ji = 1, jpi 332 ! z ind0 = 1 if sm_i le s_i_0 and 0 otherwise333 z ind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i(ji,jj,jl) ) )334 ! z ind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws335 z ind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i(ji,jj,jl) ) )336 ! If 2.sm_i GE sss_m then zindbal= 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 337 332 ! this is to force a constant salinity profile in the Baltic Sea 338 zindbal= MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )339 zalpha(ji,jj,jl) = z ind0 + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )340 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zindbal)341 END DO 342 END DO 343 END DO 344 345 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 346 341 DO jl = 1, jpl 347 342 DO jk = 1, nlay_i … … 349 344 DO ji = 1, jpi 350 345 ! ! linear profile with 0 at the surface 351 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 352 347 ! ! weighting the profile 353 348 s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 354 END DO ! ji 355 END DO ! jj 356 END DO ! jk 357 END DO ! jl 358 ! 359 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 360 357 361 358 !------------------------------------------------------- … … 363 360 !------------------------------------------------------- 364 361 365 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) 366 363 ! 367 364 sm_i(:,:,:) = 2.30_wp 368 365 ! 369 366 DO jl = 1, jpl 370 !CDIR NOVERRCHK371 367 DO jk = 1, nlay_i 372 zargtemp = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)368 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 373 369 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 374 370 s_i(:,:,jk,jl) = zsal … … 376 372 END DO 377 373 ! 378 ENDIF ! n um_sal374 ENDIF ! nn_icesal 379 375 ! 380 376 CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha ) … … 390 386 !!------------------------------------------------------------------ 391 387 INTEGER :: ji, jj, jk, jl ! dummy loop indices 392 REAL(wp) :: zindb ! - -393 388 !!------------------------------------------------------------------ 394 389 395 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 396 396 tm_i(:,:) = 0._wp 397 397 DO jl = 1, jpl … … 399 399 DO jj = 1, jpj 400 400 DO ji = 1, jpi 401 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 402 tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 403 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 404 END DO 405 END DO 406 END DO 407 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 408 409 409 410 END SUBROUTINE lim_var_icetm … … 421 422 !!------------------------------------------------------------------ 422 423 INTEGER :: ji, jj, jk, jl ! dummy loop indices 423 REAL(wp) :: zbvi, zinda, zindb ! local scalars 424 !!------------------------------------------------------------------ 425 ! 424 REAL(wp) :: zbvi ! local scalars 425 !!------------------------------------------------------------------ 426 ! 427 vt_i (:,:) = 0._wp 428 DO jl = 1, jpl 429 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 430 END DO 431 426 432 bv_i(:,:) = 0._wp 427 433 DO jl = 1, jpl … … 429 435 DO jj = 1, jpj 430 436 DO ji = 1, jpi 431 zinda = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) ) )432 z indb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) )433 zbvi = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 ) &434 & * v_i(ji,jj,jl) / REAL(nlay_i,wp)435 bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi / MAX( vt_i(ji,jj) , epsi10 )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 ) 436 442 END DO 437 443 END DO … … 452 458 ! 453 459 INTEGER :: ji, jk ! dummy loop indices 454 INTEGER :: ii, ij ! local integers455 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal ! local scalars456 REAL(wp) :: zalpha, z ind0, zind01, zindbal, 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 ! - - 457 463 ! 458 464 REAL(wp), POINTER, DIMENSION(:) :: z_slope_s 465 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 466 REAL(wp), PARAMETER :: zsi1 = 4.5_wp 459 467 !!--------------------------------------------------------------------- 460 468 … … 464 472 ! Vertically constant, constant in time 465 473 !--------------------------------------- 466 IF( n um_sal == 1 ) s_i_b(:,:) = bulk_sal474 IF( nn_icesal == 1 ) s_i_1d(:,:) = rn_icesal 467 475 468 476 !------------------------------------------------------ … … 470 478 !------------------------------------------------------ 471 479 472 IF( n um_sal == 2 ) THEN480 IF( nn_icesal == 2 ) THEN 473 481 ! 474 482 DO ji = kideb, kiut ! Slope of the linear profile zs_zero 475 z_slope_s(ji) = 2._wp * sm_i_b(ji) / MAX( epsi10 , ht_i_b(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) ) 476 485 END DO 477 486 478 487 ! Weighting factor between zs_zero and zs_inf 479 488 !--------------------------------------------- 480 dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 ) 481 dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 ) 482 dummy_fac2 = 1._wp / REAL(nlay_i,wp) 483 484 !CDIR NOVERRCHK 489 zfac0 = 1._wp / ( zsi0 - zsi1 ) 490 zfac1 = zsi1 / ( zsi1 - zsi0 ) 485 491 DO jk = 1, nlay_i 486 !CDIR NOVERRCHK487 492 DO ji = kideb, kiut 488 493 ii = MOD( npb(ji) - 1 , jpi ) + 1 489 494 ij = ( npb(ji) - 1 ) / jpi + 1 490 ! z ind0 = 1 if sm_i le s_i_0 and 0 otherwise491 z ind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_b(ji) ) )492 ! z ind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws493 z ind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) )494 ! if 2.sm_i GE sss_m then zindbal= 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 495 500 ! this is to force a constant salinity profile in the Baltic Sea 496 zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) )501 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 497 502 ! 498 zalpha = ( z ind0 + zind01 * ( sm_i_b(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zindbal)503 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 ) ) * ( 1._wp - rswitch ) 499 504 ! 500 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_ b(ji) * dummy_fac2505 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 501 506 ! weighting the profile 502 s_i_b(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_b(ji) 503 END DO ! ji 504 END DO ! jk 505 506 ENDIF ! num_sal 507 s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 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 507 514 508 515 !------------------------------------------------------- … … 510 517 !------------------------------------------------------- 511 518 512 IF( num_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 513 ! 514 sm_i_b(:) = 2.30_wp 515 ! 516 !CDIR NOVERRCHK 519 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 520 ! 521 sm_i_1d(:) = 2.30_wp 522 ! 517 523 DO jk = 1, nlay_i 518 zargtemp = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)519 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 ) ) ) ) 520 526 DO ji = kideb, kiut 521 s_i_ b(ji,jk) = zsal527 s_i_1d(ji,jk) = zsal 522 528 END DO 523 529 END DO … … 528 534 ! 529 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 530 784 531 785 #else … … 546 800 SUBROUTINE lim_var_salprof1d ! Emtpy routines 547 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 548 806 #endif 549 807
Note: See TracChangeset
for help on using the changeset viewer.
