Changeset 5034 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
- Timestamp:
- 2015-01-15T14:48:42+01: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
r4333 r5034 66 66 PUBLIC lim_var_salprof1d ! 67 67 68 REAL(wp) :: epsi10 = 1.e-10_wp ! - -69 REAL(wp) :: zzero = 0.e0 ! - -70 REAL(wp) :: zone = 1.e0 ! - -71 72 68 !!---------------------------------------------------------------------- 73 69 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 94 90 ! 95 91 INTEGER :: ji, jj, jk, jl ! dummy loop indices 96 REAL(wp) :: zinda, zindb97 92 !!------------------------------------------------------------------ 98 93 … … 113 108 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 114 109 ! 115 zinda = MAX( zzero , SIGN( zone, at_i(ji,jj) - epsi10 ) )116 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda! ice thickness110 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 111 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch ! ice thickness 117 112 END DO 118 113 END DO … … 134 129 DO jj = 1, jpj 135 130 DO ji = 1, jpi 136 zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi10 ) )137 zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) )138 131 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 139 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda ! ice salinity 140 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) * zindb ! ice age 132 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 133 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch ! ice salinity 134 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 135 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) * rswitch ! ice age 141 136 END DO 142 137 END DO … … 163 158 INTEGER :: ji, jj, jk, jl ! dummy loop indices 164 159 REAL(wp) :: zq_i, zaaa, zbbb, zccc, zdiscrim ! local scalars 165 REAL(wp) :: ztmelts, z indb, zq_s, zfac1, zfac2 ! - -160 REAL(wp) :: ztmelts, zq_s, zfac1, zfac2 ! - - 166 161 !!------------------------------------------------------------------ 167 162 … … 172 167 DO jj = 1, jpj 173 168 DO ji = 1, jpi 174 zindb= 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes175 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb176 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb177 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb169 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 170 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 171 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 172 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 178 173 END DO 179 174 END DO … … 184 179 DO jj = 1, jpj 185 180 DO ji = 1, jpi 186 zindb= 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes187 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * zindb181 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 188 183 END DO 189 184 END DO … … 205 200 DO ji = 1, jpi 206 201 ! ! Energy of melting q(S,T) [J.m-3] 207 zq_i = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)208 z indb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) ! zindb = 0 if no ice and 1 if yes209 zq_i = zq_i * unit_fac * zindb!convert units202 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 210 205 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt ! Ice layer melt temperature 211 206 ! … … 214 209 zccc = lfus * (ztmelts-rtt) 215 210 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 216 t_i(ji,jj,jk,jl) = rtt + zindb*( - zbbb - zdiscrim ) / ( 2.0 *zaaa )211 t_i(ji,jj,jk,jl) = rtt + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 217 212 t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt 218 213 END DO … … 231 226 DO ji = 1, jpi 232 227 !Energy of melting q(S,T) [J.m-3] 233 zq_s = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)234 z indb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) ) ! zindb = 0 if no ice and 1 if yes235 zq_s = zq_s * unit_fac * zindb! convert units228 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 236 231 ! 237 t_s(ji,jj,jk,jl) = rtt + zindb* ( - zfac1 * zq_s + zfac2 )232 t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 ) 238 233 t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt 239 234 END DO … … 250 245 DO jj = 1, jpj 251 246 DO ji = 1, jpi 252 zindb= ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) )253 tm_i(ji,jj) = tm_i(ji,jj) + zindb* t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) &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) & 254 249 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 255 250 END DO … … 297 292 INTEGER :: ji, jj, jk, jl ! dummy loop index 298 293 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac, zsal ! local scalar 299 REAL(wp) :: z ind0, zind01, zindbal, zargtemp , zs_zero ! - -294 REAL(wp) :: zswi0, zswi01, zswibal, zargtemp , zs_zero ! - - 300 295 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha ! 3D pointer 301 296 !!------------------------------------------------------------------ … … 320 315 DO jj = 1, jpj 321 316 DO ji = 1, jpi 322 z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( 0.01, ht_i(ji,jj,jl) )317 z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) ) 323 318 END DO 324 319 END DO … … 332 327 DO jj = 1, jpj 333 328 DO ji = 1, jpi 334 ! z ind0 = 1 if sm_i le s_i_0 and 0 otherwise335 z ind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i(ji,jj,jl) ) )336 ! z ind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws337 z ind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i(ji,jj,jl) ) )338 ! If 2.sm_i GE sss_m then z indbal = 1329 ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise 330 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 othws 332 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 = 1 339 334 ! this is to force a constant salinity profile in the Baltic Sea 340 z indbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )341 zalpha(ji,jj,jl) = z ind0 + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )342 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - z indbal )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 ) 343 338 END DO 344 339 END DO … … 392 387 !!------------------------------------------------------------------ 393 388 INTEGER :: ji, jj, jk, jl ! dummy loop indices 394 REAL(wp) :: zindb ! - -395 389 !!------------------------------------------------------------------ 396 390 … … 401 395 DO jj = 1, jpj 402 396 DO ji = 1, jpi 403 zindb= ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) )404 tm_i(ji,jj) = tm_i(ji,jj) + zindb* t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) &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) & 405 399 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 406 400 END DO … … 423 417 !!------------------------------------------------------------------ 424 418 INTEGER :: ji, jj, jk, jl ! dummy loop indices 425 REAL(wp) :: zbvi , zinda, zindb! local scalars419 REAL(wp) :: zbvi ! local scalars 426 420 !!------------------------------------------------------------------ 427 421 ! … … 431 425 DO jj = 1, jpj 432 426 DO ji = 1, jpi 433 zinda = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) ) ) 434 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 435 zbvi = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 ) & 427 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) ) ) 428 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 ) & 436 429 & * v_i(ji,jj,jl) / REAL(nlay_i,wp) 437 bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi / MAX( vt_i(ji,jj) , epsi10 ) 430 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 431 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi10 ) 438 432 END DO 439 433 END DO … … 456 450 INTEGER :: ii, ij ! local integers 457 451 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal ! local scalars 458 REAL(wp) :: zalpha, z ind0, zind01, zindbal, zs_zero ! - -452 REAL(wp) :: zalpha, zswi0, zswi01, zswibal, zs_zero ! - - 459 453 ! 460 454 REAL(wp), POINTER, DIMENSION(:) :: z_slope_s … … 466 460 ! Vertically constant, constant in time 467 461 !--------------------------------------- 468 IF( num_sal == 1 ) s_i_ b(:,:) = bulk_sal462 IF( num_sal == 1 ) s_i_1d(:,:) = bulk_sal 469 463 470 464 !------------------------------------------------------ … … 475 469 ! 476 470 DO ji = kideb, kiut ! Slope of the linear profile zs_zero 477 z_slope_s(ji) = 2._wp * sm_i_ b(ji) / MAX( 0.01 , ht_i_b(ji) )471 z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) ) 478 472 END DO 479 473 … … 490 484 ii = MOD( npb(ji) - 1 , jpi ) + 1 491 485 ij = ( npb(ji) - 1 ) / jpi + 1 492 ! z ind0 = 1 if sm_i le s_i_0 and 0 otherwise493 z ind0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_b(ji) ) )494 ! z ind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws495 z ind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) )496 ! if 2.sm_i GE sss_m then z indbal = 1486 ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise 487 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 othws 489 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 = 1 497 491 ! this is to force a constant salinity profile in the Baltic Sea 498 z indbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) )492 zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 499 493 ! 500 zalpha = ( z ind0 + zind01 * ( sm_i_b(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zindbal )494 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zswibal ) 501 495 ! 502 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_ b(ji) * dummy_fac2496 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2 503 497 ! weighting the profile 504 s_i_ b(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_b(ji)498 s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 505 499 END DO ! ji 506 500 END DO ! jk … … 514 508 IF( num_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 515 509 ! 516 sm_i_ b(:) = 2.30_wp510 sm_i_1d(:) = 2.30_wp 517 511 ! 518 512 !CDIR NOVERRCHK … … 521 515 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 522 516 DO ji = kideb, kiut 523 s_i_ b(ji,jk) = zsal517 s_i_1d(ji,jk) = zsal 524 518 END DO 525 519 END DO
Note: See TracChangeset
for help on using the changeset viewer.