- Timestamp:
- 2015-01-20T15:26:13+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4333 r5038 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 … … 35 36 36 37 PUBLIC lim_trp ! called by ice_step 37 38 REAL(wp) :: epsi10 = 1.e-10_wp39 REAL(wp) :: rzero = 0._wp40 REAL(wp) :: rone = 1._wp41 38 42 39 !! * Substitution … … 63 60 INTEGER, INTENT(in) :: kt ! number of iteration 64 61 ! 65 INTEGER :: ji, jj, jk, jl, layer! dummy loop indices62 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 66 63 INTEGER :: initad ! number of sub-timestep for the advection 67 INTEGER :: ierr ! error status 68 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 69 REAL(wp) :: zusvosn, zusvoic, zbigval ! - - 70 REAL(wp) :: zcfl , zusnit ! - - 71 REAL(wp) :: ze , zsal , zage ! - - 64 REAL(wp) :: zcfl , zusnit ! - - 72 65 ! 73 66 REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow 74 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 75 68 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 ! mass and salt flux (clem) 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold ! old ice volume... 80 ! correct ice thickness (clem) 81 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 82 REAL(wp) :: zdv, zda, zvi, zvs, zsmv 69 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume... 70 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 71 REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies 72 REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 73 ! 74 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 83 75 !!--------------------------------------------------------------------- 84 76 IF( nn_timing == 1 ) CALL timing_start('limtrp') 85 77 86 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )78 CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 87 79 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 88 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e ) 89 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 ! ------------------------------- 80 CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 81 82 CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold ) ! clem 103 83 104 84 IF( numit == nstart .AND. lwp ) THEN … … 115 95 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! 116 96 ! !-------------------------------------! 97 98 ! conservation test 99 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 100 117 101 ! mass and salt flux init (clem) 118 102 zviold(:,:,:) = v_i(:,:,:) 103 zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 104 zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 119 105 120 106 !--- Thickness correction init. (clem) ------------------------------- … … 167 153 ! ENDIF 168 154 !!gm end 169 initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )155 initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) ) 170 156 zusnit = 1.0 / REAL( initad ) 171 157 IF( zcfl > 0.5 .AND. lwp ) & … … 174 160 175 161 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 176 DO j k= 1,initad177 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area162 DO jn = 1,initad 163 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 178 164 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 179 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &165 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:), & 180 166 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 181 167 DO jl = 1, jpl 182 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---168 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 183 169 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 184 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &170 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 185 171 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 186 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---172 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 187 173 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 188 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &174 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 189 175 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 190 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---176 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 191 177 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 192 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &178 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 193 179 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 194 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---180 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 195 181 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &182 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 197 183 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 198 CALL lim_adv_x( zusnit, u_ice, rone, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---184 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 199 185 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 200 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &186 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 201 187 & 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 ---188 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 203 189 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 204 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &190 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 205 191 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 206 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 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &209 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )210 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &211 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &212 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )192 DO jk = 1, nlay_i !--- ice heat contents --- 193 CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 194 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 195 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 197 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 198 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 213 199 END DO 214 200 END DO 215 201 END DO 216 202 ELSE 217 DO j k= 1, initad218 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area203 DO jn = 1, initad 204 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 219 205 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 220 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &206 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:), & 221 207 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 222 208 DO jl = 1, jpl 223 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---209 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 224 210 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 225 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &211 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 226 212 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 227 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---213 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 228 214 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 229 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &215 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 230 216 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 231 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---217 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 232 218 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 233 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &219 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 234 220 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 235 221 236 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---222 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 237 223 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 238 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &224 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 239 225 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 240 CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---226 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 241 227 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 242 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &228 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 243 229 & 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 ---230 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 245 231 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 246 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &232 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 247 233 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 248 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 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &251 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )252 CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &253 & sxxe(:,:, layer,jl), sye (:,:,layer,jl), &254 & syye(:,:, layer,jl), sxye(:,:,layer,jl) )234 DO jk = 1, nlay_i !--- ice heat contents --- 235 CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 236 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 237 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 238 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), & 239 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 240 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 255 241 END DO 256 242 END DO … … 268 254 zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:) 269 255 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 256 ! 274 257 END DO 275 258 … … 289 272 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 290 273 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)274 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji ,jj) ) ) ) & 275 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 276 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji,jj ) ) ) ) & 277 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 295 278 END DO 296 279 END DO … … 305 288 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 306 289 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)290 pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji ,jj,jl) ) ) ) & 291 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 292 pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji,jj ,jl) ) ) ) & 293 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 311 294 END DO 312 295 END DO … … 334 317 DO jj = 1, jpj 335 318 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) )319 zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) ) 320 zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) ) 321 zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) ) 322 zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) ) 323 zs0a (ji,jj,jl) = MAX( 0._wp, zs0a (ji,jj,jl) ) 324 zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) ) 342 325 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 343 326 END DO … … 346 329 347 330 !--------------------------------------------------------- 348 ! 5.2) Snow thickness, Ice thickness, Ice concentrations331 ! 5.2) Update and mask variables 349 332 !--------------------------------------------------------- 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 333 DO jl = 1, jpl 359 334 DO jj = 1, jpj 360 335 DO ji = 1, jpi 361 zvi = zs0ice(ji,jj,jl) 362 zvs = zs0sn(ji,jj,jl) 336 rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 337 338 zvi = zs0ice(ji,jj,jl) 339 zvs = zs0sn (ji,jj,jl) 340 zes = zs0c0 (ji,jj,jl) 341 zsmv = zs0sm (ji,jj,jl) 363 342 ! 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 381 END DO 382 END DO 383 END DO 384 385 !--- Thickness correction in case too high (clem) -------------------------------------------------------- 386 CALL lim_var_glo2eqv 387 DO jl = 1, jpl 388 DO jj = 1, jpj 389 DO ji = 1, jpi 390 391 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) 395 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 396 397 zindh = 1._wp 398 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 399 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 400 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), epsi10 ) 403 ELSE 404 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), epsi10 ) 407 ENDIF 408 409 ! 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 413 ! 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 416 417 ENDIF 418 419 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 ! --- 426 DO jj = 1, jpj 427 DO ji = 1, jpi 428 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless?? 429 END DO 430 END DO 431 432 !---------------------- 433 ! 5.3) Ice properties 434 !---------------------- 435 436 zbigval = 1.e+13 437 438 DO jl = 1, jpl 439 DO jj = 1, jpj 440 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 343 ! Remove very small areas 344 v_s(ji,jj,jl) = rswitch * zs0sn (ji,jj,jl) 345 v_i(ji,jj,jl) = rswitch * zs0ice(ji,jj,jl) 346 a_i(ji,jj,jl) = rswitch * zs0a (ji,jj,jl) 347 e_s(ji,jj,1,jl) = rswitch * zs0c0 (ji,jj,jl) 450 348 ! 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 349 IF( num_sal == 2 ) THEN 453 350 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 351 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) 352 oa_i(ji,jj,jl) = MAX( rswitch * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl) 353 354 ! Update fluxes 355 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 356 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 464 357 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 358 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 359 END DO 360 END DO 361 END DO 468 362 469 363 DO jl = 1, jpl … … 471 365 DO jj = 1, jpj 472 366 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 367 rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 368 zei = zs0e(ji,jj,jk,jl) 369 e_i(ji,jj,jk,jl) = rswitch * MAX( 0._wp, zs0e(ji,jj,jk,jl) ) 370 ! Update fluxes 371 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 477 372 END DO !ji 478 373 END DO ! jj … … 480 375 END DO ! jl 481 376 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 377 !--- Thickness correction in case too high (clem) -------------------------------------------------------- 378 CALL lim_var_glo2eqv 379 DO jl = 1, jpl 380 DO jj = 1, jpj 381 DO ji = 1, jpi 382 383 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 384 zvi = v_i (ji,jj,jl) 385 zvs = v_s (ji,jj,jl) 386 zsmv = smv_i(ji,jj,jl) 387 zes = e_s (ji,jj,1,jl) 388 zei = SUM( e_i(ji,jj,1:nlay_i,jl) ) 389 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 390 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 391 392 rswitch = 1._wp 393 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 394 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 395 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 396 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 397 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 398 ELSE 399 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 400 rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 401 a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 402 ENDIF 403 404 ! small correction due to *rswitch for a_i 405 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl) 406 v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl) 407 smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl) 408 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 409 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 410 411 ! Update mass fluxes 412 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 413 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 414 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 415 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 416 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 417 ENDIF 418 END DO 419 END DO 420 END DO 421 ! ------------------------------------------------- 422 423 ! --- diags --- 489 424 DO jj = 1, jpj 490 425 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 426 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 427 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 428 429 diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 430 diag_trp_vs(ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 431 END DO 432 END DO 433 434 ! --- agglomerate variables ----------------- 435 vt_i (:,:) = 0._wp 436 vt_s (:,:) = 0._wp 437 at_i (:,:) = 0._wp 438 ! 439 DO jl = 1, jpl 440 DO jj = 1, jpj 441 DO ji = 1, jpi 442 ! 443 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 444 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 445 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 446 END DO 447 END DO 448 END DO 449 ! ------------------------------------------------- 450 451 ! open water 452 DO jj = 1, jpj 453 DO ji = 1, jpi 454 ! open water = 1 if at_i=0 455 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 456 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * zs0ow(ji,jj) 457 END DO 458 END DO 459 460 ! conservation test 461 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 504 462 505 463 ENDIF … … 536 494 END DO 537 495 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 496 ! 566 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )497 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 567 498 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 568 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )569 570 CALL wrk_dealloc( jpi, jpj,jpl,zaiold, zhimax ) ! clem499 CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 500 501 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax ) ! clem 571 502 ! 572 503 IF( nn_timing == 1 ) CALL timing_stop('limtrp')
Note: See TracChangeset
for help on using the changeset viewer.
