Changeset 4688 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
- Timestamp:
- 2014-06-25T01:39:59+02:00 (10 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4333 r4688 30 30 USE limvar ! clem for ice thickness correction 31 31 USE timing ! Timing 32 USE limcons ! conservation tests 32 33 33 34 IMPLICIT NONE … … 37 38 38 39 REAL(wp) :: epsi10 = 1.e-10_wp 39 REAL(wp) :: rzero = 0._wp 40 REAL(wp) :: rone = 1._wp 40 REAL(wp) :: epsi20 = 1.e-20_wp 41 41 42 42 !! * Substitution … … 67 67 INTEGER :: ierr ! error status 68 68 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 69 REAL(wp) :: zusvosn, zusvoic, zbigval ! - -70 69 REAL(wp) :: zcfl , zusnit ! - - 71 REAL(wp) :: z e , zsal , zage ! - -70 REAL(wp) :: zsal , zage ! - - 72 71 ! 73 72 REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow 74 73 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 75 74 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 76 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)77 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax, zchk_umax ! Check errors (C Rousset)78 75 ! mass and salt flux (clem) 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold ! old ice volume... 80 ! correct ice thickness (clem) 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 81 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 82 REAL(wp) :: zdv, zda, zvi, zvs, zsmv 78 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 79 REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 80 ! 81 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 83 82 !!--------------------------------------------------------------------- 84 83 IF( nn_timing == 1 ) CALL timing_start('limtrp') 85 84 86 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )85 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 87 86 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 88 87 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e ) 89 88 90 CALL wrk_alloc( jpi,jpj,jpl,zviold ) ! clem 91 CALL wrk_alloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem 92 93 ! ------------------------------- 94 !- check conservation (C Rousset) 95 IF( ln_limdiahsb ) THEN 96 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 97 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 98 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 99 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 100 ENDIF 101 !- check conservation (C Rousset) 102 ! ------------------------------- 89 CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold ) ! clem 103 90 104 91 IF( numit == nstart .AND. lwp ) THEN … … 115 102 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! 116 103 ! !-------------------------------------! 104 105 ! conservation test 106 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 107 117 108 ! mass and salt flux init (clem) 118 109 zviold(:,:,:) = v_i(:,:,:) 110 zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 111 zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 119 112 120 113 !--- Thickness correction init. (clem) ------------------------------- … … 167 160 ! ENDIF 168 161 !!gm end 169 initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )162 initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) ) 170 163 zusnit = 1.0 / REAL( initad ) 171 164 IF( zcfl > 0.5 .AND. lwp ) & … … 175 168 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 176 169 DO jk = 1,initad 177 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area170 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 178 171 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 179 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &172 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:), & 180 173 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 181 174 DO jl = 1, jpl 182 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---175 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 183 176 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 184 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &177 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 185 178 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 186 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---179 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 187 180 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 188 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &181 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 189 182 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 190 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---183 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 191 184 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 192 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &185 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 193 186 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 194 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---187 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 195 188 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &189 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 197 190 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 198 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---191 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 199 192 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 200 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &193 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 201 194 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 202 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---195 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 203 196 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 204 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &197 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 205 198 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 206 199 DO layer = 1, nlay_i !--- ice heat contents --- 207 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &200 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 208 201 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 209 202 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 210 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &203 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 211 204 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 212 205 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) … … 216 209 ELSE 217 210 DO jk = 1, initad 218 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area211 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 219 212 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 220 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &213 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:), & 221 214 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 222 215 DO jl = 1, jpl 223 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---216 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 224 217 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 225 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &218 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 226 219 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 227 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---220 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 228 221 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 229 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &222 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 230 223 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 231 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---224 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 232 225 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 233 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &226 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 234 227 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 235 228 236 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---229 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 237 230 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 238 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &231 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 239 232 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 240 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---233 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 241 234 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 242 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &235 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 243 236 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 244 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---237 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 245 238 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 246 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &239 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 247 240 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 248 241 DO layer = 1, nlay_i !--- ice heat contents --- 249 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &242 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 250 243 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 251 244 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 252 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &245 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 253 246 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 254 247 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) … … 268 261 zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:) 269 262 zs0a (:,:,jl) = zs0a (:,:,jl) / area(:,:) 270 zs0c0 (:,:,jl) = zs0c0 (:,:,jl) / area(:,:) 271 DO jk = 1, nlay_i 272 zs0e(:,:,jk,jl) = zs0e(:,:,jk,jl) / area(:,:) 273 END DO 263 ! 274 264 END DO 275 265 … … 289 279 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 290 280 DO ji = 1 , fs_jpim1 ! vector opt. 291 pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji ,jj) ) ) ) &292 & * ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)293 pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji,jj ) ) ) ) &294 & * ( 1._wp - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)281 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji ,jj) ) ) ) & 282 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 283 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji,jj ) ) ) ) & 284 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 295 285 END DO 296 286 END DO … … 305 295 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 306 296 DO ji = 1 , fs_jpim1 ! vector opt. 307 pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji ,jj,jl) ) ) ) &308 & * ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)309 pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji,jj ,jl) ) ) ) &310 & * ( 1._wp - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)297 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji ,jj,jl) ) ) ) & 298 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 299 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji,jj ,jl) ) ) ) & 300 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 311 301 END DO 312 302 END DO … … 334 324 DO jj = 1, jpj 335 325 DO ji = 1, jpi 336 zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) )337 zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) )338 zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) )339 zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) )340 zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl) )341 zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) )326 zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) ) 327 zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) ) 328 zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) ) 329 zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) ) 330 zs0a (ji,jj,jl) = MAX( 0._wp, zs0a (ji,jj,jl) ) 331 zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) ) 342 332 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 343 333 END DO … … 346 336 347 337 !--------------------------------------------------------- 348 ! 5.2) Snow thickness, Ice thickness, Ice concentrations338 ! 5.2) Update and mask variables 349 339 !--------------------------------------------------------- 350 DO jj = 1, jpj 351 DO ji = 1, jpi 352 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - epsi10) ) 353 zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp ) 354 ato_i(ji,jj) = zs0ow(ji,jj) 355 END DO 356 END DO 357 358 DO jl = 1, jpl ! Remove very small areas 340 DO jl = 1, jpl 359 341 DO jj = 1, jpj 360 342 DO ji = 1, jpi 361 zvi = zs0ice(ji,jj,jl) 362 zvs = zs0sn(ji,jj,jl) 343 zindb= MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 344 345 zvi = zs0ice(ji,jj,jl) 346 zvs = zs0sn (ji,jj,jl) 347 zes = zs0c0 (ji,jj,jl) 348 zsmv = zs0sm (ji,jj,jl) 363 349 ! 364 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - epsi10) ) 365 ! 366 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 367 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 368 ! 369 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 370 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 371 zindb = MAX( zindsn, zindic ) 372 ! 373 zs0a(ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration 374 a_i (ji,jj,jl) = zs0a(ji,jj,jl) 375 v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl) 376 v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl) 377 ! 378 ! Update mass fluxes (clem) 379 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 380 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 350 ! Remove very small areas 351 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 352 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 353 a_i(ji,jj,jl) = zindb * zs0a (ji,jj,jl) 354 e_s(ji,jj,1,jl) = zindb * zs0c0 (ji,jj,jl) 355 ! Ice salinity and age 356 IF( num_sal == 2 ) THEN 357 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 358 ENDIF 359 oa_i(ji,jj,jl) = MAX( zindb * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl) 360 361 ! Update fluxes 362 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 363 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 364 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 365 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 381 366 END DO 382 367 END DO 383 368 END DO 369 370 DO jl = 1, jpl 371 DO jk = 1, nlay_i 372 DO jj = 1, jpj 373 DO ji = 1, jpi 374 zindb = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 375 zei = zs0e(ji,jj,jk,jl) 376 e_i(ji,jj,jk,jl) = zindb * MAX( 0._wp, zs0e(ji,jj,jk,jl) ) 377 ! Update fluxes 378 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 379 END DO !ji 380 END DO ! jj 381 END DO ! jk 382 END DO ! jl 384 383 385 384 !--- Thickness correction in case too high (clem) -------------------------------------------------------- … … 390 389 391 390 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 392 zvi = v_i(ji,jj,jl) 393 zvs = v_s(ji,jj,jl) 394 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 391 zvi = v_i (ji,jj,jl) 392 zvs = v_s (ji,jj,jl) 393 zsmv = smv_i(ji,jj,jl) 394 zes = e_s (ji,jj,1,jl) 395 zei = SUM( e_i(ji,jj,:,jl) ) 396 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 395 397 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 396 398 … … 399 401 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 400 402 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 401 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )402 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi 10 )403 zindh = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 404 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 403 405 ELSE 404 406 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 405 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )406 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi 10 )407 zindh = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 408 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 407 409 ENDIF 408 410 409 411 ! small correction due to *zindh for a_i 410 v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) 411 v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl) 412 v_i (ji,jj,jl) = zindh * v_i (ji,jj,jl) 413 v_s (ji,jj,jl) = zindh * v_s (ji,jj,jl) 414 smv_i(ji,jj,jl) = zindh * smv_i(ji,jj,jl) 415 e_s(ji,jj,1,jl) = zindh * e_s(ji,jj,1,jl) 416 e_i(ji,jj,:,jl) = zindh * e_i(ji,jj,:,jl) 412 417 413 418 ! Update mass fluxes 414 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 415 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 419 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 420 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 421 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 422 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 423 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,:,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 416 424 417 425 ENDIF 418 426 419 427 diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 420 421 END DO 422 END DO 423 END DO 424 425 ! --- 428 diag_trp_vs(ji,jj) = diag_trp_vs(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * r1_rdtice 429 430 END DO 431 END DO 432 END DO 433 ! ------------------------------------------------- 434 435 ! --- diags --- 426 436 DO jj = 1, jpj 427 437 DO ji = 1, jpi 428 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless??429 END DO430 END DO431 432 !---------------------- 433 ! 5.3) Ice properties434 !----------------------435 436 zbigval = 1.e+13437 438 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 439 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 440 END DO 441 END DO 442 443 ! --- agglomerate variables (clem) ----------------- 444 vt_i (:,:) = 0._wp 445 vt_s (:,:) = 0._wp 446 at_i (:,:) = 0._wp 447 ! 438 448 DO jl = 1, jpl 439 449 DO jj = 1, jpj 440 450 DO ji = 1, jpi 441 zsmv = zs0sm(ji,jj,jl) 442 443 ! Switches and dummy variables 444 zusvosn = 1.0/MAX( v_s(ji,jj,jl) , epsi10 ) 445 zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi10 ) 446 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 447 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 448 zindb = MAX( zindsn, zindic ) 449 450 ! Ice salinity and age 451 !clem zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj), zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl) 452 IF( num_sal == 2 ) THEN 453 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 454 ENDIF 455 456 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) ), 0._wp ) * a_i(ji,jj,jl) 457 oa_i (ji,jj,jl) = zindic * zage 458 459 ! Snow heat content 460 ze = MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval ) 461 e_s(ji,jj,1,jl) = zindsn * ze 462 463 ! Update salt fluxes (clem) 464 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 465 END DO !ji 466 END DO !jj 467 END DO ! jl 468 469 DO jl = 1, jpl 470 DO jk = 1, nlay_i 471 DO jj = 1, jpj 472 DO ji = 1, jpi 473 ! Ice heat content 474 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 475 ze = MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval ) 476 e_i(ji,jj,jk,jl) = zindic * ze 477 END DO !ji 478 END DO ! jj 479 END DO ! jk 480 END DO ! jl 481 482 483 ! --- agglomerate variables (clem) ----------------- 484 vt_i (:,:) = 0._wp 485 vt_s (:,:) = 0._wp 486 at_i (:,:) = 0._wp 487 ! 488 DO jl = 1, jpl 451 ! 452 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 453 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 454 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 455 END DO 456 END DO 457 END DO 458 ! ------------------------------------------------- 459 460 ! open water 489 461 DO jj = 1, jpj 490 462 DO ji = 1, jpi 491 ! 492 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 493 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 494 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 495 ! 496 zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi10 ) ) 497 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda ! ice thickness 498 END DO 499 END DO 500 END DO 501 ! ------------------------------------------------- 502 503 463 ! open water = 1 if at_i=0 464 zindb = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 465 ato_i(ji,jj) = zindb + (1._wp - zindb ) * zs0ow(ji,jj) 466 END DO 467 END DO 468 469 ! conservation test 470 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 504 471 505 472 ENDIF … … 536 503 END DO 537 504 ENDIF 538 ! -------------------------------539 !- check conservation (C Rousset)540 IF( ln_limdiahsb ) THEN541 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b542 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b543 544 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice545 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )546 547 zchk_vmin = glob_min(v_i)548 zchk_amax = glob_max(SUM(a_i,dim=3))549 zchk_amin = glob_min(a_i)550 zchk_umax = glob_max(SQRT(u_ice**2 + v_ice**2))551 552 IF(lwp) THEN553 IF ( ABS( zchk_v_i ) > 1.e-5 ) THEN554 WRITE(numout,*) 'violation volume [m3/day] (limtrp) = ',(zchk_v_i * rday)555 WRITE(numout,*) 'u_ice max [m/s] (limtrp) = ',zchk_umax556 WRITE(numout,*) 'number of time steps (limtrp) =',kt557 ENDIF558 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday)559 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limtrp) = ',(zchk_vmin * 1.e-3)560 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limtrp) = ',zchk_amin561 ENDIF562 ENDIF563 !- check conservation (C Rousset)564 ! -------------------------------565 505 ! 566 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )506 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 567 507 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 568 508 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e ) 569 509 570 CALL wrk_dealloc( jpi, jpj,jpl,zaiold, zhimax ) ! clem510 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax ) ! clem 571 511 ! 572 512 IF( nn_timing == 1 ) CALL timing_stop('limtrp')
Note: See TracChangeset
for help on using the changeset viewer.
