Changeset 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
- Timestamp:
- 2015-04-13T15:08:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r4688 r5208 66 66 PUBLIC lim_var_salprof1d ! 67 67 68 REAL(wp) :: epsi10 = 1.e-10_wp ! - -69 70 68 !!---------------------------------------------------------------------- 71 69 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 92 90 ! 93 91 INTEGER :: ji, jj, jk, jl ! dummy loop indices 94 REAL(wp) :: zinda, zindb95 92 !!------------------------------------------------------------------ 96 93 … … 111 108 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 112 109 ! 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 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 115 112 END DO 116 113 END DO … … 132 129 DO jj = 1, jpj 133 130 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 131 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 137 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda ! ice salinity 138 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 139 136 END DO 140 137 END DO … … 161 158 INTEGER :: ji, jj, jk, jl ! dummy loop indices 162 159 REAL(wp) :: zq_i, zaaa, zbbb, zccc, zdiscrim ! local scalars 163 REAL(wp) :: ztmelts, z indb, zq_s, zfac1, zfac2 ! - -160 REAL(wp) :: ztmelts, zq_s, zfac1, zfac2 ! - - 164 161 !!------------------------------------------------------------------ 165 162 … … 170 167 DO jj = 1, jpj 171 168 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) , epsi10 ) * zindb174 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb175 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 176 173 END DO 177 174 END DO … … 182 179 DO jj = 1, jpj 183 180 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 yes185 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 186 183 END DO 187 184 END DO … … 203 200 DO ji = 1, jpi 204 201 ! ! 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 yes206 zq_i = zindb* e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)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) 207 204 zq_i = zq_i * unit_fac !convert units 208 205 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt ! Ice layer melt temperature … … 212 209 zccc = lfus * (ztmelts-rtt) 213 210 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 214 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 ) 215 212 t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt 216 213 END DO … … 229 226 DO ji = 1, jpi 230 227 !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 yes232 zq_s = zindb* e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)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) 233 230 zq_s = zq_s * unit_fac ! convert units 234 231 ! 235 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 ) 236 233 t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt 237 234 END DO … … 248 245 DO jj = 1, jpj 249 246 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) &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) & 252 249 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 253 250 END DO … … 295 292 INTEGER :: ji, jj, jk, jl ! dummy loop index 296 293 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac, zsal ! local scalar 297 REAL(wp) :: z ind0, zind01, zindbal, zargtemp , zs_zero ! - -294 REAL(wp) :: zswi0, zswi01, zswibal, zargtemp , zs_zero ! - - 298 295 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha ! 3D pointer 299 296 !!------------------------------------------------------------------ … … 330 327 DO jj = 1, jpj 331 328 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 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 337 334 ! this is to force a constant salinity profile in the Baltic Sea 338 z indbal = 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 - 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 ) 341 338 END DO 342 339 END DO … … 390 387 !!------------------------------------------------------------------ 391 388 INTEGER :: ji, jj, jk, jl ! dummy loop indices 392 REAL(wp) :: zindb ! - -393 389 !!------------------------------------------------------------------ 394 390 … … 399 395 DO jj = 1, jpj 400 396 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) &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) & 403 399 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ) ) 404 400 END DO … … 421 417 !!------------------------------------------------------------------ 422 418 INTEGER :: ji, jj, jk, jl ! dummy loop indices 423 REAL(wp) :: zbvi , zinda, zindb! local scalars419 REAL(wp) :: zbvi ! local scalars 424 420 !!------------------------------------------------------------------ 425 421 ! … … 429 425 DO jj = 1, jpj 430 426 DO ji = 1, jpi 431 zinda = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) ) ) 432 zindb = ( 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 ) & 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 ) & 434 429 & * 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 ) 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 ) 436 432 END DO 437 433 END DO … … 454 450 INTEGER :: ii, ij ! local integers 455 451 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal ! local scalars 456 REAL(wp) :: zalpha, z ind0, zind01, zindbal, zs_zero ! - -452 REAL(wp) :: zalpha, zswi0, zswi01, zswibal, zs_zero ! - - 457 453 ! 458 454 REAL(wp), POINTER, DIMENSION(:) :: z_slope_s … … 464 460 ! Vertically constant, constant in time 465 461 !--------------------------------------- 466 IF( num_sal == 1 ) s_i_ b(:,:) = bulk_sal462 IF( num_sal == 1 ) s_i_1d(:,:) = bulk_sal 467 463 468 464 !------------------------------------------------------ … … 473 469 ! 474 470 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) )471 z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) ) 476 472 END DO 477 473 … … 488 484 ii = MOD( npb(ji) - 1 , jpi ) + 1 489 485 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 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 495 491 ! this is to force a constant salinity profile in the Baltic Sea 496 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) ) ) 497 493 ! 498 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 ) 499 495 ! 500 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 501 497 ! weighting the profile 502 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) 503 499 END DO ! ji 504 500 END DO ! jk … … 512 508 IF( num_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 513 509 ! 514 sm_i_ b(:) = 2.30_wp510 sm_i_1d(:) = 2.30_wp 515 511 ! 516 512 !CDIR NOVERRCHK … … 519 515 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 520 516 DO ji = kideb, kiut 521 s_i_ b(ji,jk) = zsal517 s_i_1d(ji,jk) = zsal 522 518 END DO 523 519 END DO
Note: See TracChangeset
for help on using the changeset viewer.