- Timestamp:
- 2014-05-12T22:46:18+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4332 r4634 37 37 38 38 REAL(wp) :: epsi10 = 1.e-10_wp 39 REAL(wp) :: rzero = 0._wp 40 REAL(wp) :: rone = 1._wp 39 REAL(wp) :: epsi20 = 1.e-20_wp 41 40 42 41 !! * Substitution … … 67 66 INTEGER :: ierr ! error status 68 67 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 69 REAL(wp) :: zusvosn, zusvoic, zbigval ! - -70 68 REAL(wp) :: zcfl , zusnit ! - - 71 REAL(wp) :: z e , zsal , zage ! - -69 REAL(wp) :: zsal , zage ! - - 72 70 ! 73 71 REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow 74 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 75 73 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)74 REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset) 77 75 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax, zchk_umax ! Check errors (C Rousset) 78 76 ! mass and salt flux (clem) 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold ! old ice volume... 80 ! correct ice thickness (clem) 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 81 78 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 82 REAL(wp) :: zdv, zda, zvi, zvs, zsmv 79 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 80 REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 83 81 !!--------------------------------------------------------------------- 84 82 IF( nn_timing == 1 ) CALL timing_start('limtrp') 85 83 86 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )84 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 87 85 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 88 86 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e ) 89 87 90 CALL wrk_alloc( jpi,jpj,jpl,zviold ) ! clem 91 CALL wrk_alloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem 88 CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold ) ! clem 92 89 93 90 ! ------------------------------- 94 91 !- check conservation (C Rousset) 95 92 IF( ln_limdiahsb ) THEN 96 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:) , dim=3 ) * area(:,:) * tms(:,:) )93 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 97 94 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(:,:) ) 95 zchk_e_i_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 96 zchk_fw_b = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 97 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 98 zchk_ft_b = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 100 99 ENDIF 101 100 !- check conservation (C Rousset) … … 117 116 ! mass and salt flux init (clem) 118 117 zviold(:,:,:) = v_i(:,:,:) 118 zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 119 zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 119 120 120 121 !--- Thickness correction init. (clem) ------------------------------- … … 167 168 ! ENDIF 168 169 !!gm end 169 initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )170 initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) ) 170 171 zusnit = 1.0 / REAL( initad ) 171 172 IF( zcfl > 0.5 .AND. lwp ) & … … 175 176 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 176 177 DO jk = 1,initad 177 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area178 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 178 179 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 179 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &180 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:), & 180 181 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 181 182 DO jl = 1, jpl 182 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---183 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 183 184 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 184 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &185 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 185 186 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 186 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---187 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 187 188 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 188 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &189 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 189 190 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 190 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---191 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 191 192 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 192 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &193 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 193 194 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 194 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---195 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 195 196 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &197 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 197 198 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 198 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---199 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 199 200 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 200 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &201 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 201 202 & 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 ---203 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 203 204 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 204 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &205 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 205 206 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 206 207 DO layer = 1, nlay_i !--- ice heat contents --- 207 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &208 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 208 209 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 209 210 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 210 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &211 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 211 212 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 212 213 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) … … 216 217 ELSE 217 218 DO jk = 1, initad 218 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area219 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 219 220 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 220 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &221 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:), & 221 222 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 222 223 DO jl = 1, jpl 223 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---224 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 224 225 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 225 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &226 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 226 227 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 227 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---228 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 228 229 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 229 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &230 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 230 231 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 231 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---232 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 232 233 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 233 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &234 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 234 235 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 235 236 236 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---237 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 237 238 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 238 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &239 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 239 240 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 240 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---241 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 241 242 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 242 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &243 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 243 244 & 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 ---245 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 245 246 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 246 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &247 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 247 248 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 248 249 DO layer = 1, nlay_i !--- ice heat contents --- 249 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &250 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 250 251 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 251 252 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 252 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &253 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 253 254 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 254 255 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) … … 268 269 zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:) 269 270 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 271 ! 274 272 END DO 275 273 … … 289 287 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 290 288 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)289 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji ,jj) ) ) ) & 290 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 291 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji,jj ) ) ) ) & 292 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 295 293 END DO 296 294 END DO … … 305 303 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 306 304 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)305 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji ,jj,jl) ) ) ) & 306 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 307 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji,jj ,jl) ) ) ) & 308 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 311 309 END DO 312 310 END DO … … 334 332 DO jj = 1, jpj 335 333 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) )334 zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) ) 335 zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) ) 336 zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) ) 337 zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) ) 338 zs0a (ji,jj,jl) = MAX( 0._wp, zs0a (ji,jj,jl) ) 339 zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) ) 342 340 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 343 341 END DO … … 346 344 347 345 !--------------------------------------------------------- 348 ! 5.2) Snow thickness, Ice thickness, Ice concentrations346 ! 5.2) Update and mask variables 349 347 !--------------------------------------------------------- 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 348 DO jl = 1, jpl 359 349 DO jj = 1, jpj 360 350 DO ji = 1, jpi 361 zvi = zs0ice(ji,jj,jl) 362 zvs = zs0sn(ji,jj,jl) 351 zindb= MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 352 353 zvi = zs0ice(ji,jj,jl) 354 zvs = zs0sn (ji,jj,jl) 355 zes = zs0c0 (ji,jj,jl) 356 zsmv = zs0sm (ji,jj,jl) 363 357 ! 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 358 ! Remove very small areas 359 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 360 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 361 a_i(ji,jj,jl) = zindb * zs0a (ji,jj,jl) 362 e_s(ji,jj,1,jl) = zindb * zs0c0 (ji,jj,jl) 363 ! Ice salinity and age 364 IF( num_sal == 2 ) THEN 365 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 366 ENDIF 367 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) 368 369 ! Update fluxes 370 wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 371 wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 372 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 373 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 374 END DO 382 375 END DO 383 376 END DO 377 378 DO jl = 1, jpl 379 DO jk = 1, nlay_i 380 DO jj = 1, jpj 381 DO ji = 1, jpi 382 zindb = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 383 zei = zs0e(ji,jj,jk,jl) 384 e_i(ji,jj,jk,jl) = zindb * MAX( 0._wp, zs0e(ji,jj,jk,jl) ) 385 ! Update fluxes 386 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 387 END DO !ji 388 END DO ! jj 389 END DO ! jk 390 END DO ! jl 384 391 385 392 !--- Thickness correction in case too high (clem) -------------------------------------------------------- … … 390 397 391 398 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) 399 zvi = v_i (ji,jj,jl) 400 zvs = v_s (ji,jj,jl) 401 zsmv = smv_i(ji,jj,jl) 402 zes = e_s (ji,jj,1,jl) 403 zei = SUM( e_i(ji,jj,:,jl) ) 404 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 395 405 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 396 406 … … 399 409 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 400 410 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 )411 zindh = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 412 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 403 413 ELSE 404 414 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 )415 zindh = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 416 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 407 417 ENDIF 408 418 409 419 ! 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) 420 v_i (ji,jj,jl) = zindh * v_i (ji,jj,jl) 421 v_s (ji,jj,jl) = zindh * v_s (ji,jj,jl) 422 smv_i(ji,jj,jl) = zindh * smv_i(ji,jj,jl) 423 e_s(ji,jj,1,jl) = zindh * e_s(ji,jj,1,jl) 424 e_i(ji,jj,:,jl) = zindh * e_i(ji,jj,:,jl) 412 425 413 426 ! 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 427 wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 428 wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 429 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 430 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 431 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 432 417 433 ENDIF 418 434 419 435 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 ! --- 436 diag_trp_vs(ji,jj) = diag_trp_vs(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * r1_rdtice 437 438 END DO 439 END DO 440 END DO 441 ! ------------------------------------------------- 442 443 ! --- diags --- 426 444 DO jj = 1, jpj 427 445 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 446 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 447 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 448 END DO 449 END DO 450 451 ! --- agglomerate variables (clem) ----------------- 452 vt_i (:,:) = 0._wp 453 vt_s (:,:) = 0._wp 454 at_i (:,:) = 0._wp 455 ! 438 456 DO jl = 1, jpl 439 457 DO jj = 1, jpj 440 458 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 459 ! 460 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 461 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 462 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 463 END DO 464 END DO 465 END DO 466 ! ------------------------------------------------- 467 468 ! open water 489 469 DO jj = 1, jpj 490 470 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 471 ! open water = 1 if at_i=0 472 zindb = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 473 ato_i(ji,jj) = zindb + (1._wp - zindb ) * zs0ow(ji,jj) 474 END DO 475 END DO 504 476 505 477 ENDIF … … 539 511 !- check conservation (C Rousset) 540 512 IF( ln_limdiahsb ) THEN 541 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 542 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 513 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 514 zchk_fw = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 515 zchk_ft = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 543 516 544 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 545 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 517 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 518 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 519 zchk_e_i = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 546 520 547 521 zchk_vmin = glob_min(v_i) … … 551 525 552 526 IF(lwp) THEN 553 IF ( ABS( zchk_v_i ) > 1.e- 5) THEN554 WRITE(numout,*) 'violation volume [ m3/day] (limtrp) = ',(zchk_v_i * rday)527 IF ( ABS( zchk_v_i ) > 1.e-4 ) THEN 528 WRITE(numout,*) 'violation volume [kg/day] (limtrp) = ',(zchk_v_i * rday) 555 529 WRITE(numout,*) 'u_ice max [m/s] (limtrp) = ',zchk_umax 556 530 WRITE(numout,*) 'number of time steps (limtrp) =',kt 557 531 ENDIF 558 532 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday) 533 IF ( ABS( zchk_e_i ) > 1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J] (limtrp) = ',(zchk_e_i) 559 534 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limtrp) = ',(zchk_vmin * 1.e-3) 560 535 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limtrp) = ',zchk_amin … … 564 539 ! ------------------------------- 565 540 ! 566 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )541 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 567 542 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 568 543 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e ) 569 544 570 CALL wrk_dealloc( jpi, jpj,jpl,zaiold, zhimax ) ! clem545 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax ) ! clem 571 546 ! 572 547 IF( nn_timing == 1 ) CALL timing_stop('limtrp')
Note: See TracChangeset
for help on using the changeset viewer.
