Changeset 5078
- Timestamp:
- 2015-02-11T16:15:11+01:00 (9 years ago)
- Location:
- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
- Files:
-
- 13 edited
-
LIM_SRC_3/ice.F90 (modified) (3 diffs)
-
LIM_SRC_3/limcons.F90 (modified) (1 diff)
-
LIM_SRC_3/limistate.F90 (modified) (3 diffs)
-
LIM_SRC_3/limitd_me.F90 (modified) (1 diff)
-
LIM_SRC_3/limthd.F90 (modified) (3 diffs)
-
LIM_SRC_3/limthd_dh.F90 (modified) (11 diffs)
-
LIM_SRC_3/limthd_dif.F90 (modified) (8 diffs)
-
LIM_SRC_3/limthd_ent.F90 (modified) (1 diff)
-
LIM_SRC_3/limthd_lac.F90 (modified) (4 diffs)
-
LIM_SRC_3/limvar.F90 (modified) (13 diffs)
-
LIM_SRC_3/limwri.F90 (modified) (1 diff)
-
OPA_SRC/DOM/phycst.F90 (modified) (2 diffs)
-
OPA_SRC/SBC/sbcice_lim.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r5076 r5078 165 165 166 166 ! !!** ice-thickness distribution namelist (namiceitd) ** 167 INTEGER , PUBLIC :: jpl !: number of ice categories168 INTEGER , PUBLIC :: nlay_i !: number of ice layers169 INTEGER , PUBLIC :: nlay_s !: number of snow layers170 167 INTEGER , PUBLIC :: nn_catbnd !: categories distribution following: tanh function (1), or h^(-alpha) function (2) 171 168 REAL(wp), PUBLIC :: rn_himean !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only) … … 223 220 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( rn_ecc * rn_ecc ) 224 221 REAL(wp), PUBLIC :: rhoco !: = rau0 * cio 225 222 REAL(wp), PUBLIC :: r1_nlay_i !: 1 / nlay_i 223 REAL(wp), PUBLIC :: r1_nlay_s !: 1 / nlay_s 224 ! 226 225 ! !!** switch for presence of ice or not 227 226 REAL(wp), PUBLIC :: rswitch 228 227 ! 229 228 ! !!** define some parameters 230 229 REAL(wp), PUBLIC, PARAMETER :: epsi06 = 1.e-06_wp !: small number … … 368 367 !!-------------------------------------------------------------------------- 369 368 ! !!: ** Namelist namicerun read in sbc_lim_init ** 370 CHARACTER(len=32) , PUBLIC :: cn_icerst_in !: suffix of ice restart name (input) 371 CHARACTER(len=32) , PUBLIC :: cn_icerst_out !: suffix of ice restart name (output) 372 LOGICAL , PUBLIC :: ln_limdyn !: flag for ice dynamics (T) or not (F) 373 LOGICAL , PUBLIC :: ln_nicep !: flag for sea-ice points output (T) or not (F) 374 REAL(wp) , PUBLIC :: rn_amax !: maximum ice concentration 369 INTEGER , PUBLIC :: jpl !: number of ice categories 370 INTEGER , PUBLIC :: nlay_i !: number of ice layers 371 INTEGER , PUBLIC :: nlay_s !: number of snow layers 372 CHARACTER(len=32), PUBLIC :: cn_icerst_in !: suffix of ice restart name (input) 373 CHARACTER(len=32), PUBLIC :: cn_icerst_out !: suffix of ice restart name (output) 374 LOGICAL , PUBLIC :: ln_limdyn !: flag for ice dynamics (T) or not (F) 375 LOGICAL , PUBLIC :: ln_nicep !: flag for sea-ice points output (T) or not (F) 376 REAL(wp) , PUBLIC :: rn_amax !: maximum ice concentration 375 377 ! 376 378 !!-------------------------------------------------------------------------- -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r5076 r5078 210 210 & * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw 211 211 212 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs /rhoic )212 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs * r1_rhoic ) 213 213 214 214 zei = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + & -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r5067 r5078 326 326 ! recompute ht_i, ht_s avoiding out of bounds values 327 327 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 328 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic /rhosn )328 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 329 329 330 330 ! ice volume, salt content, age content … … 347 347 348 348 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 349 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) /nlay_s349 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 350 350 END DO ! ji 351 351 END DO ! jj … … 368 368 369 369 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 370 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) /nlay_i370 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 371 371 END DO ! ji 372 372 END DO ! jj -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r5070 r5078 1124 1124 ! clem: if sst>0, then ersw <0 (is that possible?) 1125 1125 zsstK = sst_m(ji,jj) + rt0 1126 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i )1126 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) * r1_nlay_i 1127 1127 1128 1128 ! heat flux to the ocean -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5076 r5078 337 337 DO jl = 1, jpl 338 338 DO jk = 1, nlay_i 339 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) / REAL( nlay_i )339 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i 340 340 END DO 341 341 END DO … … 347 347 DO jl = 1, jpl 348 348 DO jk = 1, nlay_s 349 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) / REAL( nlay_s )349 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s 350 350 END DO 351 351 END DO … … 479 479 ! Conversion q(S,T) -> T (second order equation) 480 480 zaaa = cpic 481 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) /rhoic - lfus481 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) * r1_rhoic - lfus 482 482 zccc = lfus * ( ztmelts - rtt ) 483 483 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5070 r5078 145 145 DO jk = 1, nlay_i 146 146 DO ji = kideb, kiut 147 h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i )147 h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 148 148 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 149 149 ENDDO … … 188 188 ! 189 189 DO ji = kideb, kiut 190 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )190 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 191 191 END DO 192 192 ! … … 199 199 DO jk = 1, nlay_i 200 200 DO ji = kideb, kiut 201 zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i )201 zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 202 202 zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 203 203 END DO … … 228 228 ! thickness change 229 229 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji) 230 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice /rhosn230 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 231 231 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 232 232 zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) … … 255 255 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 256 256 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 257 zh_s (ji) = ht_s_1d(ji) / REAL( nlay_s )257 zh_s (ji) = ht_s_1d(ji) * r1_nlay_s 258 258 259 259 ENDIF … … 311 311 DO ji = kideb, kiut 312 312 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 313 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )313 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 314 314 END DO ! ji 315 315 … … 335 335 DO jk = 1, nlay_i 336 336 DO ji = kideb, kiut 337 zEi = - q_i_1d(ji,jk) / rhoic! Specific enthalpy of layer k [J/kg, <0]337 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 338 338 339 339 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer k [K] … … 345 345 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 346 346 347 zdeltah(ji,jk) = - zfmdt / rhoic! Melt of layer jk [m, <0]347 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] 348 348 349 349 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] … … 506 506 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 507 507 508 zEi = - q_i_1d(ji,jk) / rhoic! Specific enthalpy of melting ice (J/kg, <0)508 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 509 509 510 510 !!zEw = rcp * ( t_i_1d(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) … … 535 535 ELSE !!! Basal melting 536 536 537 zEi = - q_i_1d(ji,jk) /rhoic ! Specific enthalpy of melting ice (J/kg, <0)538 539 zEw = rcp * ( ztmelts - rtt ) ! Specific enthalpy of meltwater (J/kg, <0)540 541 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)542 543 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0)544 545 zdeltah(ji,jk) = - zfmdt / rhoic! Gross thickness change537 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 538 539 zEw = rcp * ( ztmelts - rtt ) ! Specific enthalpy of meltwater (J/kg, <0) 540 541 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 542 543 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 544 545 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change 546 546 547 547 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change … … 627 627 ! Salinity of snow ice 628 628 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 629 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) /rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)629 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 630 630 631 631 ! entrapment during snow ice formation -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5076 r5078 179 179 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 180 180 DO ji = kideb, kiut 181 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i )+ &182 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ))181 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 182 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 183 183 END DO 184 184 … … 189 189 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ! is there snow or not 190 190 ! layer thickness 191 zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i )192 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )191 zh_i(ji) = ht_i_1d(ji) * r1_nlay_i 192 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 193 193 END DO 194 194 … … 202 202 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 203 203 DO ji = kideb , kiut 204 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s )204 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 205 205 END DO 206 206 END DO … … 208 208 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 209 209 DO ji = kideb , kiut 210 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i )210 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 211 211 END DO 212 212 END DO … … 576 576 !------------------------------------------------------------------------------| 577 577 ! 578 IF ( t_su_1d(ji) < rtt) THEN578 IF ( t_su_1d(ji) < rtt ) THEN 579 579 ! 580 580 !------------------------------------------------------------------------------| … … 600 600 !!case of only one layer in the ice (surface & ice equations are altered) 601 601 602 IF ( nlay_i.eq.1) THEN602 IF ( nlay_i == 1 ) THEN 603 603 ztrid(ji,numeqmin(ji),1) = 0.0 604 604 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0) * 2.0 … … 631 631 632 632 !!case of only one layer in the ice (surface & ice equations are altered) 633 IF ( nlay_i.eq.1) THEN633 IF ( nlay_i == 1 ) THEN 634 634 ztrid(ji,numeqmin(ji),1) = 0.0 635 635 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) … … 760 760 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 761 761 DO ji = kideb, kiut 762 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i )+ &763 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ))762 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 763 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 764 764 IF( t_su_1d(ji) < rtt ) THEN ! case T_su < 0degC 765 765 zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r5064 r5078 102 102 ! new layer thickesses 103 103 DO ji = kideb, kiut 104 zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) / REAL( nlay_i )104 zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i 105 105 ENDDO 106 106 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r5070 r5078 362 362 DO ji = 1, nbpac 363 363 364 zEi = - ze_newice(ji) / rhoic! specific enthalpy of forming ice [J/kg]365 366 zEw = rcp * ( t_bo_1d(ji) - rt0 ) 364 zEi = - ze_newice(ji) * r1_rhoic ! specific enthalpy of forming ice [J/kg] 365 366 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 367 367 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 368 368 … … 371 371 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 372 372 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 373 zv_newice(ji) = - zfmdt /rhoic373 zv_newice(ji) = - zfmdt * r1_rhoic 374 374 375 375 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux … … 458 458 DO jk = 1, nlay_i 459 459 DO ji = 1, nbpac 460 h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i )460 h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 461 461 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 462 462 END DO … … 525 525 DO ji = 1, jpi 526 526 ! heat content in J/m2 527 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / REAL( nlay_i ,wp )527 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 528 528 END DO 529 529 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r5067 r5078 196 196 ! 197 197 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 198 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i /rhoic - lfus198 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i * r1_rhoic - lfus 199 199 zccc = lfus * (ztmelts-rtt) 200 200 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) … … 280 280 !!------------------------------------------------------------------ 281 281 INTEGER :: ji, jj, jk, jl ! dummy loop index 282 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac, zsal283 REAL(wp) :: zswi0, zswi01, z swibal, zargtemp , zs_zero282 REAL(wp) :: zfac0, zfac1, zsal 283 REAL(wp) :: zswi0, zswi01, zargtemp , zs_zero 284 284 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha 285 285 REAL(wp), PARAMETER :: zsi0 = 3.5_wp … … 311 311 END DO 312 312 ! 313 dummy_fac0 = 1._wp / ( zsi0 - zsi1 ) ! Weighting factor between zs_zero and zs_inf314 dummy_fac1 = zsi1/ ( zsi1 - zsi0 )313 zfac0 = 1._wp / ( zsi0 - zsi1 ) ! Weighting factor between zs_zero and zs_inf 314 zfac1 = zsi1 / ( zsi1 - zsi0 ) 315 315 ! 316 316 zalpha(:,:,:) = 0._wp … … 322 322 ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 323 323 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i(ji,jj,jl) ) ) 324 ! If 2.sm_i GE sss_m then zswibal= 1324 ! If 2.sm_i GE sss_m then rswitch = 1 325 325 ! this is to force a constant salinity profile in the Baltic Sea 326 zswibal= MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )327 zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )328 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal)329 END DO 330 END DO 331 END DO 332 333 dummy_fac = 1._wp / REAL( nlay_i )! Computation of the profile326 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 327 zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 ) 328 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch ) 329 END DO 330 END DO 331 END DO 332 333 ! Computation of the profile 334 334 DO jl = 1, jpl 335 335 DO jk = 1, nlay_i … … 337 337 DO ji = 1, jpi 338 338 ! ! linear profile with 0 at the surface 339 zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac339 zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i 340 340 ! ! weighting the profile 341 341 s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) … … 357 357 DO jl = 1, jpl 358 358 DO jk = 1, nlay_i 359 zargtemp = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)359 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 360 360 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 361 361 s_i(:,:,jk,jl) = zsal … … 387 387 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 388 388 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 389 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ))389 & * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 ) 390 390 END DO 391 391 END DO … … 417 417 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) ) ) 418 418 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 ) & 419 & * v_i(ji,jj,jl) / REAL(nlay_i,wp)419 & * v_i(ji,jj,jl) * r1_nlay_i 420 420 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 421 421 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi10 ) … … 439 439 INTEGER :: ji, jk ! dummy loop indices 440 440 INTEGER :: ii, ij ! local integers 441 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal ! local scalars442 REAL(wp) :: zalpha, zswi0, zswi01, zs wibal, zs_zero ! - -441 REAL(wp) :: zfac0, zfac1, zargtemp, zsal ! local scalars 442 REAL(wp) :: zalpha, zswi0, zswi01, zs_zero ! - - 443 443 ! 444 444 REAL(wp), POINTER, DIMENSION(:) :: z_slope_s … … 466 466 ! Weighting factor between zs_zero and zs_inf 467 467 !--------------------------------------------- 468 dummy_fac0 = 1._wp / ( zsi0 - zsi1 ) 469 dummy_fac1 = zsi1 / ( zsi1 - zsi0 ) 470 dummy_fac2 = 1._wp / REAL(nlay_i,wp) 471 468 zfac0 = 1._wp / ( zsi0 - zsi1 ) 469 zfac1 = zsi1 / ( zsi1 - zsi0 ) 472 470 DO jk = 1, nlay_i 473 471 DO ji = kideb, kiut … … 478 476 ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 479 477 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 480 ! if 2.sm_i GE sss_m then zswibal= 1478 ! if 2.sm_i GE sss_m then rswitch = 1 481 479 ! this is to force a constant salinity profile in the Baltic Sea 482 zswibal= MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )480 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 483 481 ! 484 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zswibal)482 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 ) ) * ( 1.0 - rswitch ) 485 483 ! 486 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2484 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 487 485 ! weighting the profile 488 486 s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) … … 501 499 ! 502 500 DO jk = 1, nlay_i 503 zargtemp = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)501 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 504 502 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 505 503 DO ji = kideb, kiut … … 734 732 ! recompute ht_i, ht_s avoiding out of bounds values 735 733 zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) 736 zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic /rhosn )734 zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn ) 737 735 ENDIF 738 736 ENDDO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r5076 r5078 266 266 zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* & 267 267 ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - epsi06 ) ) * & 268 rswitch /nlay_i268 rswitch * r1_nlay_i 269 269 END DO 270 270 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r4990 r5078 82 82 REAL(wp), PUBLIC :: xlic = 300.33e+6_wp !: volumetric latent heat fusion of ice [J/m3] 83 83 REAL(wp), PUBLIC :: xsn = 2.8e+6_wp !: volumetric latent heat of sublimation of snow [J/m3] 84 #endif 85 #if defined key_lim3 86 REAL(wp), PUBLIC :: r1_rhoic !: 1 / rhoic 87 REAL(wp), PUBLIC :: r1_rhosn !: 1 / rhosn 84 88 #endif 85 89 !!---------------------------------------------------------------------- … … 166 170 lfus = xlsn / rhosn ! latent heat of fusion of fresh ice 167 171 #endif 168 172 #if defined key_lim3 173 r1_rhoic = 1._wp / rhoic 174 r1_rhosn = 1._wp / rhosn 175 #endif 169 176 IF(lwp) THEN 170 177 WRITE(numout,*) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r5070 r5078 391 391 rdt_ice = nn_fsbc * rdttra(1) 392 392 r1_rdtice = 1._wp / rdt_ice 393 394 ! inverse of nlay_i and nlay_s 395 r1_nlay_i = 1._wp / REAL( nlay_i, wp ) 396 r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 393 397 ! 394 398 END SUBROUTINE ice_run
Note: See TracChangeset
for help on using the changeset viewer.
