- Timestamp:
- 2017-09-01T15:49:35+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8426 r8486 31 31 !! History : - ! 2006-01 (M. Vancoppenolle) Original code 32 32 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 33 !! 3.5 ! 2012 (M. Vancoppenolle) add ice_var_itd 34 !! 3.6 ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd 33 35 !!---------------------------------------------------------------------- 34 36 #if defined key_lim3 35 37 !!---------------------------------------------------------------------- 36 38 !! 'key_lim3' LIM3 sea-ice model 39 !!---------------------------------------------------------------------- 40 !! ice_var_agg : integrate variables over layers and categories 41 !! ice_var_glo2eqv : transform from VGLO to VEQV 42 !! ice_var_eqv2glo : transform from VEQV to VGLO 43 !! ice_var_salprof : salinity profile in the ice 44 !! ice_var_bv : brine volume 45 !! ice_var_salprof1d : salinity profile in the ice 1D 46 !! ice_var_zapsmall : remove very small area and volume 47 !! ice_var_itd : convert 1-cat to multiple cat 37 48 !!---------------------------------------------------------------------- 38 49 USE par_oce ! ocean parameters … … 59 70 60 71 !!---------------------------------------------------------------------- 61 !! NEMO/ LIM3 3.5 , UCL - NEMO Consortium (2011)72 !! NEMO/ICE 4.0 , NEMO Consortium (2017) 62 73 !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $ 63 74 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 69 80 !! *** ROUTINE ice_var_agg *** 70 81 !! 71 !! ** Purpose : aggregates ice-thickness-category variables to all-ice variables 72 !! i.e. it turns VGLO into VAGG 73 !! ** Method : 74 !! 75 !! ** Arguments : n = 1, at_i vt_i only 76 !! n = 2 everything 77 !! 78 !! note : you could add an argument when you need only at_i, vt_i 79 !! and when you need everything 80 !!------------------------------------------------------------------ 81 INTEGER, INTENT( in ) :: kn ! =1 at_i & vt only ; = what is needed 82 ! 83 INTEGER :: ji, jj, jk, jl ! dummy loop indices 84 !!------------------------------------------------------------------ 85 86 ! integrated values 87 vt_i (:,:) = SUM( v_i, dim=3 ) 88 vt_s (:,:) = SUM( v_s, dim=3 ) 89 at_i (:,:) = SUM( a_i, dim=3 ) 82 !! ** Purpose : aggregates ice-thickness-category variables to 83 !! all-ice variables, i.e. it turns VGLO into VAGG 84 !!------------------------------------------------------------------ 85 INTEGER, INTENT( in ) :: kn ! =1 state variables only 86 ! ! >1 state variables + others 87 ! 88 INTEGER :: ji, jj, jk, jl ! dummy loop indices 89 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i, z1_vt_i 90 !!------------------------------------------------------------------ 91 ! 92 ! ! integrated values 93 vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 ) 94 vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 ) 95 at_i(:,:) = SUM( a_i(:,:,:) , dim=3 ) 90 96 et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 91 97 et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 92 98 93 99 ! MV MP 2016 94 IF ( ln_pnd ) THEN 95 at_ip(:,:) = SUM( a_ip , dim=3 )96 vt_ip(:,:) = SUM( v_ip , dim=3 )100 IF ( ln_pnd ) THEN ! Melt pond 101 at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) 102 vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) 97 103 ENDIF 98 104 ! END MP 2016 99 105 100 ! open water fraction 101 DO jj = 1, jpj 106 DO jj = 1, jpj ! open water fraction 102 107 DO ji = 1, jpi 103 108 ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp ) 104 109 END DO 105 110 END DO 111 !!gm I think this should do the work : 112 ! ato_i(:,:) = MAX( 1._wp - at_i(:,:), 0._wp ) 113 !!gm end 106 114 107 115 IF( kn > 1 ) THEN 108 109 ! mean ice/snow thickness 110 DO jj = 1, jpj 111 DO ji = 1, jpi 112 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 113 htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 114 htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 115 ENDDO 116 ENDDO 117 118 ! mean temperature (K), salinity and age 119 smt_i(:,:) = 0._wp 120 tm_i(:,:) = 0._wp 121 tm_su(:,:) = 0._wp 122 tm_si(:,:) = 0._wp 123 om_i (:,:) = 0._wp 124 DO jl = 1, jpl 125 126 DO jj = 1, jpj 127 DO ji = 1, jpi 128 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 129 tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) 130 tm_si(ji,jj) = tm_si(ji,jj) + rswitch * ( t_si(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) 131 om_i (ji,jj) = om_i (ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 ) 132 END DO 133 END DO 134 135 DO jk = 1, nlay_i 136 DO jj = 1, jpj 137 DO ji = 1, jpi 138 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 139 tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) & 140 & / MAX( vt_i(ji,jj) , epsi10 ) 141 smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 142 & / MAX( vt_i(ji,jj) , epsi10 ) 143 END DO 144 END DO 145 END DO 146 END DO 147 tm_i = tm_i + rt0 148 tm_su = tm_su + rt0 149 tm_si = tm_si + rt0 150 ! 116 ! 117 ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) ) 118 WHERE( at_i(:,:) > epsi10 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:) 119 ELSEWHERE ; z1_at_i(:,:) = 0._wp 120 END WHERE 121 WHERE( vt_i(:,:) > epsi10 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:) 122 ELSEWHERE ; z1_vt_i(:,:) = 0._wp 123 END WHERE 124 ! 125 ! ! mean ice/snow thickness 126 htm_i(:,:) = vt_i(:,:) * z1_at_i(:,:) 127 htm_s(:,:) = vt_s(:,:) * z1_at_i(:,:) 128 ! 129 ! ! mean temperature (K), salinity and age 130 tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(ji,jj) 131 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(ji,jj) 132 om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(ji,jj) 133 ! 134 tm_i (:,:) = r1_nlay_i * SUM( SUM( t_i(:,:,:,:) * v_i(:,:,:), dim=4 ) , dim=3 ) * z1_vt_i(:,:) 135 smt_i(:,:) = r1_nlay_i * SUM( SUM( s_i(:,:,:,:) * v_i(:,:,:), dim=4 ) , dim=3 ) * z1_vt_i(:,:) 136 ! 137 !!gm QUESTION 1 : why salinity is named smt_i and not just sm_i ? since the 4D field is named s_i. (NB for temp: tm_i, t_i) 138 ! 139 DEALLOCATE( z1_at_i , z1_vt_i ) 151 140 ENDIF 152 141 ! … … 158 147 !! *** ROUTINE ice_var_glo2eqv *** 159 148 !! 160 !! ** Purpose : computes equivalent variables as function of global variables161 !! i.e. it turns VGLO into VEQV149 !! ** Purpose : computes equivalent variables as function of 150 !! global variables, i.e. it turns VGLO into VEQV 162 151 !!------------------------------------------------------------------ 163 152 INTEGER :: ji, jj, jk, jl ! dummy loop indices 164 REAL(wp) :: ze_i, zaaa, zbbb, zccc, zdiscrim ! local scalars 165 REAL(wp) :: ztmelts, ze_s, zfac1, zfac2 ! - - 166 !!------------------------------------------------------------------ 153 REAL(wp) :: ze_i, z1_CpR, zdiscrim, zaaa, z1_2aaa ! local scalars 154 REAL(wp) :: ze_s, zL_Cp , ztmelts , zbbb, zccc ! - - 155 REAL(wp) :: zhmax, z1_zhmax, zsm_i, zcpMcpic, zt_i ! - - 156 REAL(wp) :: zlay_i, zlay_s ! - - 157 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i 158 !!------------------------------------------------------------------ 159 160 !!gm Question 1: here use epsi20 , in ice_var_agg it is epsi10 which is used... why ?? 161 162 !!gm Question 2: It is possible to define existence of sea-ice in a common way between 163 !! ice area and ice volume ? 164 !! the idea is to be able to define one for all at the begining of this routine 165 !! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 ) 167 166 168 167 !------------------------------------------------------- 169 168 ! Ice thickness, snow thickness, ice salinity, ice age 170 169 !------------------------------------------------------- 171 DO jl = 1, jpl 172 DO jj = 1, jpj 173 DO ji = 1, jpi 174 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes 175 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 176 END DO 177 END DO 178 END DO 179 ! Force the upper limit of ht_i to always be < hi_max (99 m). 180 DO jj = 1, jpj 181 DO ji = 1, jpi 182 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) ) 183 ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) ) 184 a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch 185 END DO 186 END DO 187 188 DO jl = 1, jpl 189 DO jj = 1, jpj 190 DO ji = 1, jpi 191 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes 192 ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 193 o_i(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 194 END DO 195 END DO 196 END DO 170 ! !--- inverse of the ice area 171 WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 172 ELSEWHERE ; z1_a_i(:,:,:) = 0._wp 173 END WHERE 174 ! 175 ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness 176 177 zhmax = hi_max(jpl) 178 z1_zhmax = 1._wp / hi_max(jpl) 179 WHERE( ht_i(:,:,jpl) > zhmax ) !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area 180 ht_i (:,:,jpl) = zhmax 181 a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 182 z1_a_i(:,:,jpl) = zhmax / v_i(:,:,jpl) ! NB: v_i always /=0 as ht_i > hi_max 183 END WHERE 184 185 ht_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) !--- snow thickness 197 186 198 IF( nn_icesal == 2 )THEN 199 DO jl = 1, jpl 200 DO jj = 1, jpj 201 DO ji = 1, jpi 202 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice and 1 if yes 203 sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch 204 ! ! bounding salinity 205 sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin ) 206 END DO 207 END DO 208 END DO 187 o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) !--- ice age 188 189 IF( nn_icesal == 2 )THEN !--- salinity (with a minimum value imposed everywhere) 190 WHERE( v_i(:,:,:) > epsi20 ) ; sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) ) 191 ELSEWHERE ; sm_i(:,:,:) = rn_simin 192 END WHERE 209 193 ENDIF 210 194 … … 212 196 213 197 !------------------- 214 ! Ice temperature s198 ! Ice temperature [K] (with a minimum value (rt0 - 100.) imposed everywhere) 215 199 !------------------- 200 zlay_i = REAL( nlay_i , wp ) ! number of layers 201 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 202 z1_2aaa = 1._wp / 2._wp *zaaa 203 zcpMcpic = rcp - cpic 216 204 DO jl = 1, jpl 217 205 DO jk = 1, nlay_i 218 206 DO jj = 1, jpj 219 207 DO ji = 1, jpi 220 ! ! Energy of melting e(S,T) [J.m-3] 221 rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 222 ze_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 223 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0 ! Ice layer melt temperature 224 ! 225 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 226 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + ze_i * r1_rhoic - lfus 227 zccc = lfus * (ztmelts-rt0) 228 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 229 t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 230 t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) ) ! -100 < t_i < ztmelts 208 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 209 ! 210 ze_i = e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 211 ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C] 212 ! 213 zbbb = zcpMcpic * ztmelts + ze_i * r1_rhoic - lfus 214 zccc = lfus * ztmelts 215 zdiscrim = SQRT( MAX( zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 216 zt_i = - ( zbbb + zdiscrim ) * z1_2aaa ! [C] 217 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( zt_i , ztmelts ) ) + rt0 ! [K] with bounds: -100 < zt_i < ztmelts 218 ! 219 ELSE !--- no ice 220 t_i(ji,jj,jk,jl) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) 221 !!clem: I think we should impose rt0 instead 222 ENDIF 231 223 END DO 232 224 END DO … … 235 227 236 228 !-------------------- 237 ! Snow temperature s229 ! Snow temperature [K] (with a minimum value (rt0 - 100.) imposed everywhere) 238 230 !-------------------- 239 zfac1 = 1._wp / ( rhosn * cpic ) 240 zfac2 = lfus / cpic 241 DO jl = 1, jpl 242 DO jk = 1, nlay_s 243 DO jj = 1, jpj 244 DO ji = 1, jpi 245 !Energy of melting q(S,T) [J.m-3] 246 rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 247 ze_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 248 ! 249 t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * ze_s + zfac2 ) 250 t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) ) ! -100 < t_s < rt0 251 END DO 252 END DO 253 END DO 231 z1_CpR = 1._wp / ( cpic * rhosn ) 232 zL_Cp = lfus / cpic 233 zlay_s = REAL( nlay_s , wp ) 234 DO jk = 1, nlay_s 235 WHERE( v_s(:,:,:) > epsi20 ) !--- icy area 236 t_s(:,:,jk,:) = MAX( -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp ) ) + rt0 237 ELSEWHERE !--- no ice 238 t_s(:,:,jk,:) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) 239 END WHERE 254 240 END DO 241 !!gm perhaps better like this (?) : (if this compile, yes! one test instead of nlay_s tests) 242 ! WHERE( v_s(:,:,:) > epsi20 ) !--- icy area 243 ! DO jk = 1, nlay_s 244 ! t_s(:,:,jk,:) = MAX( -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp ) ) + rt0 245 ! END DO 246 ! ELSEWHERE !--- no ice 247 ! DO jk = 1, nlay_s 248 ! t_s(:,:,jk,:) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) 249 ! END DO 250 ! END WHERE 251 !!gm end better ? 255 252 256 253 ! integrated values … … 259 256 at_i (:,:) = SUM( a_i, dim=3 ) 260 257 261 262 258 ! MV MP 2016 259 ! probably should resum for melt ponds ??? 263 260 264 261 ! … … 270 267 !! *** ROUTINE ice_var_eqv2glo *** 271 268 !! 272 !! ** Purpose : computes global variables as function of equivalent variables 273 !! i.e. it turns VEQV into VGLO 274 !! ** Method : 275 !! 276 !! ** History : (01-2006) Martin Vancoppenolle, UCL-ASTR 277 !!------------------------------------------------------------------ 278 ! 279 v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:) 280 v_s(:,:,:) = ht_s(:,:,:) * a_i(:,:,:) 269 !! ** Purpose : computes global variables as function of 270 !! equivalent variables, i.e. it turns VEQV into VGLO 271 !!------------------------------------------------------------------ 272 ! 273 v_i (:,:,:) = ht_i(:,:,:) * a_i(:,:,:) 274 v_s (:,:,:) = ht_s(:,:,:) * a_i(:,:,:) 281 275 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 282 276 ! … … 300 294 !!------------------------------------------------------------------ 301 295 INTEGER :: ji, jj, jk, jl ! dummy loop index 302 REAL(wp) :: z fac0, zfac1, zsal303 REAL(wp) :: z swi0, zswi01, zargtemp , zs_zero304 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_slope_s, zalpha296 REAL(wp) :: zsal, z1_dS 297 REAL(wp) :: zargtemp , zs0, zs 298 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only 305 299 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 306 300 REAL(wp), PARAMETER :: zsi1 = 4.5_wp 307 301 !!------------------------------------------------------------------ 308 302 309 !--------------------------------------- 310 ! Vertically constant, constant in time 311 !--------------------------------------- 312 IF( nn_icesal == 1 ) THEN 303 !!gm much much more secure to defined when reading nn_icesal in the namelist integers =1, 2, 3 with explicit names 304 !! for example np_Scst_noProfile = 1 ; np_Svar_linProfile = 2 ; np_Scst_fixProfile 305 306 !!gm Question: Remove the option 3 ? How many years since it last use ? 307 308 SELECT CASE ( nn_icesal ) 309 ! 310 ! !---------------------------------------! 311 CASE( 1 ) ! constant salinity in time and space ! 312 ! !---------------------------------------! 313 313 s_i (:,:,:,:) = rn_icesal 314 314 sm_i(:,:,:) = rn_icesal 315 ENDIF316 317 !-----------------------------------318 ! Salinity profile, varying in time319 !-----------------------------------320 IF( nn_icesal == 2 ) THEN315 ! 316 ! !---------------------------------------------! 317 CASE( 2 ) ! time varying salinity with linear profile ! 318 ! !---------------------------------------------! 319 ! 320 ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) ) 321 321 ! 322 322 DO jk = 1, nlay_i 323 323 s_i(:,:,jk,:) = sm_i(:,:,:) 324 324 END DO 325 ! 326 DO jl = 1, jpl ! Slope of the linear profile 327 DO jj = 1, jpj 328 DO ji = 1, jpi 329 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) ) 330 z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) ) 331 END DO 332 END DO 333 END DO 334 ! 335 zfac0 = 1._wp / ( zsi0 - zsi1 ) ! Weighting factor between zs_zero and zs_inf 336 zfac1 = zsi1 / ( zsi1 - zsi0 ) 337 ! 338 zalpha(:,:,:) = 0._wp 325 ! ! Slope of the linear profile 326 WHERE( ht_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * sm_i(:,:,:) / ht_i(:,:,:) 327 ELSEWHERE ; z_slope_s(:,:,:) = 0._wp 328 END WHERE 329 ! 330 z1_dS = 1._wp / ( zsi1 - zsi0 ) 339 331 DO jl = 1, jpl 340 332 DO jj = 1, jpj 341 333 DO ji = 1, jpi 342 ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise 343 zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i(ji,jj,jl) ) ) 344 ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 345 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i(ji,jj,jl) ) ) 346 ! If 2.sm_i GE sss_m then rswitch = 1 347 ! this is to force a constant salinity profile in the Baltic Sea 348 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 349 zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 ) 350 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch ) 334 zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - sm_i(ji,jj,jl) ) * z1_dS , 1._wp ) ) 335 ! ! force a constant profile when SSS too low (Baltic Sea) 336 IF( 2._wp * sm_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp 351 337 END DO 352 338 END DO … … 358 344 DO jj = 1, jpj 359 345 DO ji = 1, jpi 360 ! ! linear profile with 0 at the surface 361 zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i 362 ! ! weighting the profile 363 s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 364 ! ! bounding salinity 365 s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) ) 346 ! ! linear profile with 0 surface value 347 zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i 348 zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) ! weighting the profile 349 s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) 366 350 END DO 367 351 END DO … … 369 353 END DO 370 354 ! 371 ENDIF ! nn_icesal 372 373 !------------------------------------------------------- 374 ! Vertically varying salinity profile, constant in time 375 !------------------------------------------------------- 376 377 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 355 DEALLOCATE( z_slope_s , zalpha ) 356 ! 357 ! !-------------------------------------------! 358 CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile 359 ! !-------------------------------------------! (mean = 2.30) 378 360 ! 379 361 sm_i(:,:,:) = 2.30_wp 362 !!gm Remark: if we keep the case 3, then compute an store one for all time-step 363 !! a array S_prof(1:nlay_i) containing the calculation and just do: 364 ! DO jk = 1, nlay_i 365 ! s_i(:,:,jk,:) = S_prof(jk) 366 ! END DO 367 !!gm end 380 368 ! 381 369 DO jl = 1, jpl 382 370 DO jk = 1, nlay_i 383 371 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 384 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 385 s_i(:,:,jk,jl) = zsal 386 END DO 387 END DO 388 ! 389 ENDIF ! nn_icesal 372 s_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 373 END DO 374 END DO 375 ! 376 END SELECT 390 377 ! 391 378 END SUBROUTINE ice_var_salprof … … 405 392 !!------------------------------------------------------------------ 406 393 ! 407 bvm_i(:,:) = 0._wp 394 !!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done 395 !! instead of setting everything to zero as just below 408 396 bv_i (:,:,:) = 0._wp 409 397 DO jl = 1, jpl 410 398 DO jk = 1, nlay_i 411 DO jj = 1, jpj 412 DO ji = 1, jpi 413 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) ) 414 bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i & 415 & / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) 416 END DO 417 END DO 418 END DO 419 420 DO jj = 1, jpj 421 DO ji = 1, jpi 422 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 423 bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 ) 424 END DO 399 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) 400 bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 401 END WHERE 425 402 END DO 426 403 END DO 427 ! 404 WHERE( vt_i(:,:) > epsi20 ) bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 405 ELSEWHERE bvm_i(:,:) = 0._wp 406 END WHERE 407 ! 428 408 END SUBROUTINE ice_var_bv 429 409 … … 437 417 !!------------------------------------------------------------------- 438 418 INTEGER :: ji, jk ! dummy loop indices 439 REAL(wp) :: z fac0, zfac1, zargtemp, zsal! local scalars440 REAL(wp) :: zalpha, zs wi0, zswi01, zs_zero! - -441 ! 442 REAL(wp), DIMENSION(jpij) :: z_slope_s419 REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars 420 REAL(wp) :: zalpha, zs, zs0 ! - - 421 ! 422 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s ! 443 423 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 444 424 REAL(wp), PARAMETER :: zsi1 = 4.5_wp 445 425 !!--------------------------------------------------------------------- 446 447 !--------------------------------------- 448 ! Vertically constant, constant in time 449 !--------------------------------------- 450 IF( nn_icesal == 1 ) s_i_1d(:,:) = rn_icesal 451 452 !------------------------------------------------------ 453 ! Vertically varying salinity profile, varying in time 454 !------------------------------------------------------ 455 456 IF( nn_icesal == 2 ) THEN 457 ! 458 DO ji = 1, nidx ! Slope of the linear profile zs_zero 426 ! 427 SELECT CASE ( nn_icesal ) 428 ! 429 ! !---------------------------------------! 430 CASE( 1 ) ! constant salinity in time and space ! 431 ! !---------------------------------------! 432 s_i_1d(:,:) = rn_icesal 433 ! 434 ! !---------------------------------------------! 435 CASE( 2 ) ! time varying salinity with linear profile ! 436 ! !---------------------------------------------! 437 ! 438 ALLOCATE( z_slope_s(jpij) ) 439 ! 440 DO ji = 1, nidx ! Slope of the linear profile 459 441 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 460 442 z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) ) 461 443 END DO 462 444 463 ! Weighting factor between zs_zero and zs_inf 464 !--------------------------------------------- 465 zfac0 = 1._wp / ( zsi0 - zsi1 ) 466 zfac1 = zsi1 / ( zsi1 - zsi0 ) 445 z1_dS = 1._wp / ( zsi1 - zsi0 ) 467 446 DO jk = 1, nlay_i 468 447 DO ji = 1, nidx 469 ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise 470 zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i_1d(ji) ) ) 471 ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 472 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 473 ! if 2.sm_i GE sss_m then rswitch = 1 474 ! this is to force a constant salinity profile in the Baltic Sea 475 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_1d(ji) ) ) 448 zalpha = MAX( 0._wp , MIN( ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp ) ) 449 IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) ) zalpha = 0._wp ! cst profile when SSS too low (Baltic Sea) 476 450 ! 477 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 ) ) * ( 1._wp - rswitch ) 478 ! 479 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 480 ! weighting the profile 481 s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 482 ! bounding salinity 483 s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) ) 484 END DO 485 END DO 486 487 ENDIF 488 489 !------------------------------------------------------- 490 ! Vertically varying salinity profile, constant in time 491 !------------------------------------------------------- 492 493 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 451 zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i ! linear profile with 0 surface value 452 zs = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji) ! weighting the profile 453 s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) ! bounding salinity 454 END DO 455 END DO 456 ! 457 DEALLOCATE( z_slope_s ) 458 459 ! !-------------------------------------------! 460 CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile 461 ! !-------------------------------------------! (mean = 2.30) 494 462 ! 495 463 sm_i_1d(:) = 2.30_wp 496 464 ! 465 !!gm cf remark in ice_var_salprof routine, CASE( 3 ) 497 466 DO jk = 1, nlay_i 498 467 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i … … 503 472 END DO 504 473 ! 505 END IF474 END SELECT 506 475 ! 507 476 END SUBROUTINE ice_var_salprof1d 477 508 478 509 479 SUBROUTINE ice_var_zapsmall … … 512 482 !! 513 483 !! ** Purpose : Remove too small sea ice areas and correct fluxes 514 !!515 !! history : LIM3.5 - 01-2014 (C. Rousset) original code516 484 !!------------------------------------------------------------------- 517 485 INTEGER :: ji, jj, jl, jk ! dummy loop indices 518 486 REAL(wp) :: zsal, zvi, zvs, zei, zes, zvp 519 487 !!------------------------------------------------------------------- 520 DO jl = 1, jpl 521 488 ! 489 DO jl = 1, jpl !== loop over the categories ==! 490 ! 522 491 !----------------------------------------------------------------- 523 492 ! Zap ice energy and use ocean heat to melt ice … … 526 495 DO jj = 1 , jpj 527 496 DO ji = 1 , jpi 497 !!gm I think we can do better/faster : to be discussed 528 498 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 529 499 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch … … 588 558 589 559 ! to be sure that at_i is the sum of a_i(jl) 590 at_i (:,:) = 0._wp 591 DO jl = 1, jpl 560 at_i (:,:) = a_i(:,:,1) 561 vt_i (:,:) = v_i(:,:,1) 562 DO jl = 2, jpl 592 563 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 564 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 593 565 END DO 594 566 595 ! open water = 1 if at_i=0 567 ! open water = 1 if at_i=0 (no re-calculation of ato_i here) 596 568 DO jj = 1, jpj 597 569 DO ji = 1, jpi … … 600 572 END DO 601 573 END DO 602 603 574 ! 604 575 END SUBROUTINE ice_var_zapsmall 576 605 577 606 578 SUBROUTINE ice_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i ) … … 633 605 !! 634 606 !! (Example of application: BDY forcings when input are cell averaged) 635 !!636 607 !!------------------------------------------------------------------- 637 !! History : LIM3.5 - 2012 (M. Vancoppenolle) Original code638 !! 2014 (C. Rousset) Rewriting639 !!-------------------------------------------------------------------640 !! Local variables641 608 INTEGER :: ji, jk, jl ! dummy loop indices 642 609 INTEGER :: ijpij, i_fill, jl0 … … 645 612 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables 646 613 INTEGER , DIMENSION(4) :: itest 647 614 !!------------------------------------------------------------------- 615 648 616 !-------------------------------------------------------------------- 649 617 ! initialisation of variables 650 618 !-------------------------------------------------------------------- 651 ijpij = SIZE( zhti,1)619 ijpij = SIZE( zhti , 1 ) 652 620 zht_i(1:ijpij,1:jpl) = 0._wp 653 621 zht_s(1:ijpij,1:jpl) = 0._wp … … 766 734 zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn ) 767 735 ENDIF 768 END DO769 END DO736 END DO 737 END DO 770 738 ! 771 739 END SUBROUTINE ice_var_itd 772 773 740 774 741 #else … … 776 743 !! Default option Dummy module NO LIM3 sea-ice model 777 744 !!---------------------------------------------------------------------- 778 CONTAINS779 SUBROUTINE ice_var_itd780 END SUBROUTINE ice_var_itd781 745 #endif 782 746
Note: See TracChangeset
for help on using the changeset viewer.
