Changeset 8504 for branches/2017/dev_r8183_ICEMODEL
- Timestamp:
- 2017-09-06T17:47:28+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 3 edited
-
iceadv.F90 (modified) (12 diffs)
-
iceadv_prather.F90 (modified) (5 diffs)
-
iceadv_umx.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
r8500 r8504 36 36 PUBLIC ice_adv ! called by icestp 37 37 38 INTEGER :: ncfl ! number of ice time step with CFL>1/239 40 38 !! * Substitution 41 39 # include "vectopt_loop_substitute.h90" … … 51 49 !! *** ROUTINE ice_adv *** 52 50 !! 53 !! ** purpose : advection /diffusion processof sea ice51 !! ** purpose : advection of sea ice 54 52 !! 55 53 !! ** method : variables included in the process are scalar, … … 72 70 REAL(wp), DIMENSION(jpi,jpj) :: zatold, zeiold, zesold, zsmvold 73 71 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhimax, zviold, zvsold 74 ! --- ultimate macho only --- !75 REAL(wp) :: zdt76 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zudy, zvdx, zcu_box, zcv_box77 ! --- prather only --- !78 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zarea79 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0opw80 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0es , z0smi , z0oi81 ! MV MP 201682 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ap , z0vp83 ! END MV MP 201684 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z0ei85 72 !!--------------------------------------------------------------------- 86 73 IF( nn_timing == 1 ) CALL timing_start('iceadv') 87 74 88 75 IF( kt == nit000 .AND. lwp ) THEN 89 WRITE(numout,*)'' 90 WRITE(numout,*)'iceadv : sea-ice advection' 91 WRITE(numout,*)'~~~~~~~' 92 ncfl = 0 ! nb of time step with CFL > 1/2 76 WRITE(numout,*) 77 WRITE(numout,*) 'iceadv: sea-ice advection' 78 WRITE(numout,*) '~~~~~~' 93 79 ENDIF 94 80 95 81 CALL ice_var_agg( 1 ) ! integrated values + ato_i 96 97 !-------------------------------------!98 ! Advection of sea ice properties !99 !-------------------------------------!100 82 101 83 ! conservation test … … 103 85 104 86 ! store old values for diag 105 zviold = v_i106 zvsold = v_s107 zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 )108 zeiold (:,:) = et_i109 zesold (:,:) = et_s110 111 ! --- Thickness correction init. --- !87 zviold (:,:,:) = v_i(:,:,:) 88 zvsold (:,:,:) = v_s(:,:,:) 89 zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 ) 90 zeiold (:,:) = et_i(:,:) 91 zesold (:,:) = et_s(:,:) 92 93 ! Thickness correction init. 112 94 zatold(:,:) = at_i 113 95 WHERE( a_i(:,:,:) >= epsi20 ) … … 119 101 END WHERE 120 102 121 ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- !103 ! Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick 122 104 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 123 105 DO jl = 1, jpl … … 130 112 END DO 131 113 CALL lbc_lnk( zhimax(:,:,:), 'T', 1. ) 132 133 ! --- If ice drift field is too fast, use an appropriate time step for advection --- ! 134 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability 135 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 136 IF( lk_mpp ) CALL mpp_max( zcfl ) 137 138 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 139 ELSE ; initad = 1 ; zusnit = 1.0_wp 140 ENDIF 141 142 !! IF( zcfl > 0.5_wp .AND. lwp ) THEN 143 !! ncfl = ncfl + 1 144 !! IF( ncfl > 0 ) THEN 145 !! WRITE(cltmp,'(i6.1)') ncfl 146 !! CALL ctl_warn( 'ice_adv: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 147 !! ENDIF 148 !! ENDIF 149 114 115 !---------- 116 ! Advection 117 !---------- 150 118 SELECT CASE ( nn_limadv ) 151 152 !=============================! 153 CASE ( 0 ) !== Ultimate-MACHO scheme ==! 154 !=============================! 155 156 ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 157 158 IF( kt == nit000 .AND. lwp ) THEN 159 WRITE(numout,*)'' 160 WRITE(numout,*)'ice_adv_umx : Ultimate-MACHO advection scheme' 161 WRITE(numout,*)'~~~~~~~~~~~' 162 ENDIF 163 ! 164 zdt = rdt_ice / REAL(initad) 165 166 ! transport 167 zudy(:,:) = u_ice(:,:) * e2u(:,:) 168 zvdx(:,:) = v_ice(:,:) * e1v(:,:) 169 170 ! define velocity for advection: u*grad(H) 171 DO jj = 2, jpjm1 172 DO ji = fs_2, fs_jpim1 173 IF ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp 174 ELSEIF( u_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = u_ice(ji-1,jj) 175 ELSE ; zcu_box(ji,jj) = u_ice(ji ,jj) 176 ENDIF 177 178 IF ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp 179 ELSEIF( v_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = v_ice(ji,jj-1) 180 ELSE ; zcv_box(ji,jj) = v_ice(ji,jj ) 181 ENDIF 182 END DO 183 END DO 184 185 ! advection 186 DO jt = 1, initad 187 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) ) ! Open water area 188 DO jl = 1, jpl 189 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) ) ! Ice area 190 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) ) ! Ice volume 191 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) ) ! Salt content 192 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) ) ! Age content 193 DO jk = 1, nlay_i 194 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) ) ! Ice heat content 195 END DO 196 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) ) ! Snow volume 197 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) ) ! Snow heat content 198 ! MV MP 2016 199 IF ( nn_pnd_scheme > 0 ) THEN 200 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_ip(:,:,jl) ) ! Melt pond fraction 201 CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_ip(:,:,jl) ) ! Melt pond volume 202 ENDIF 203 ! END MV MP 2016 204 END DO 205 END DO 206 ! 207 at_i(:,:) = a_i(:,:,1) ! total ice fraction 208 DO jl = 2, jpl 209 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 210 END DO 211 ! 212 DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 213 214 !=============================! 215 CASE ( -1 ) !== Prather scheme ==! 216 !=============================! 217 218 ALLOCATE( zarea(jpi,jpj) , z0opw(jpi,jpj, 1 ) , & 219 & z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) , & 220 & z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) , & 221 & z0ei (jpi,jpj,nlay_i,jpl) ) 222 223 IF( kt == nit000 .AND. lwp ) THEN 224 WRITE(numout,*)'' 225 WRITE(numout,*)'ice_adv_xy : Prather advection scheme' 226 WRITE(numout,*)'~~~~~~~~~~~' 227 ENDIF 228 229 zarea(:,:) = e1e2t(:,:) 230 231 !------------------------- 232 ! transported fields 233 !------------------------- 234 z0opw(:,:,1) = ato_i(:,:) * e1e2t(:,:) ! Open water area 235 DO jl = 1, jpl 236 z0snw(:,:,jl) = v_s (:,:, jl) * e1e2t(:,:) ! Snow volume 237 z0ice(:,:,jl) = v_i (:,:, jl) * e1e2t(:,:) ! Ice volume 238 z0ai (:,:,jl) = a_i (:,:, jl) * e1e2t(:,:) ! Ice area 239 z0smi(:,:,jl) = smv_i(:,:, jl) * e1e2t(:,:) ! Salt content 240 z0oi (:,:,jl) = oa_i (:,:, jl) * e1e2t(:,:) ! Age content 241 z0es (:,:,jl) = e_s (:,:,1,jl) * e1e2t(:,:) ! Snow heat content 242 DO jk = 1, nlay_i 243 z0ei(:,:,jk,jl) = e_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 244 END DO 245 ! MV MP 2016 246 IF ( nn_pnd_scheme > 0 ) THEN 247 z0ap(:,:,jl) = a_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 248 z0vp(:,:,jl) = v_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 249 ENDIF 250 ! END MV MP 2016 251 END DO 252 253 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 254 DO jt = 1, initad 255 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 256 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 257 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), & 258 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 259 DO jl = 1, jpl 260 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 261 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 262 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & 263 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 264 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 265 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 266 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & 267 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 268 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 269 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 270 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & 271 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 272 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 273 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 274 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & 275 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 276 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 277 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 278 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & 279 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 280 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 281 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 282 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & 283 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 284 DO jk = 1, nlay_i !--- ice heat contents --- 285 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 286 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 287 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 288 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 289 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 290 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 291 END DO 292 ! MV MP 2016 293 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & !--- melt pond fraction -- 294 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 295 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & 296 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 297 CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & !--- melt pond volume -- 298 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 299 CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & 300 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 301 ! END MV MP 2016 302 END DO 303 END DO 304 ELSE 305 DO jt = 1, initad 306 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 307 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 308 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), & 309 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 310 DO jl = 1, jpl 311 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 312 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 313 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & 314 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 315 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 316 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 317 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & 318 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 319 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 320 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 321 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & 322 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 323 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 324 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 325 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & 326 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 327 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 328 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 329 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & 330 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 331 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 332 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 333 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & 334 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 335 DO jk = 1, nlay_i !--- ice heat contents --- 336 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 337 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 338 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 339 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 340 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 341 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 342 END DO 343 ! MV MP 2016 344 IF ( nn_pnd_scheme > 0 ) THEN 345 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & !--- melt pond fraction --- 346 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 347 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & 348 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 349 CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & !--- melt pond volume --- 350 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 351 CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & 352 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 353 ENDIF 354 ! END MV MP 2016 355 END DO 356 END DO 357 ENDIF 358 359 !------------------------------------------- 360 ! Recover the properties from their contents 361 !------------------------------------------- 362 ato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) 363 DO jl = 1, jpl 364 v_i (:,:, jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) 365 v_s (:,:, jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) 366 smv_i(:,:, jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) 367 oa_i (:,:, jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) 368 a_i (:,:, jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) 369 e_s (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:) 370 DO jk = 1, nlay_i 371 e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) 372 END DO 373 ! MV MP 2016 374 IF ( nn_pnd_scheme > 0 ) THEN 375 a_ip (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) 376 v_ip (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) 377 ENDIF 378 ! END MV MP 2016 379 END DO 380 ! 381 at_i(:,:) = a_i(:,:,1) ! total ice fraction 382 DO jl = 2, jpl 383 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 384 END DO 385 ! 386 DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei ) 387 ! 119 CASE ( 0 ) !-- Ultimate-MACHO scheme 120 CALL ice_adv_umx( kt, u_ice, v_ice, & 121 & ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 122 123 CASE ( -1 ) !-- Prather scheme 124 CALL ice_adv_prather( kt, u_ice, v_ice, & 125 & ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 126 388 127 END SELECT 389 390 ! --- diags --- 128 129 ! total ice fraction 130 at_i(:,:) = a_i(:,:,1) 131 DO jl = 2, jpl 132 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 133 END DO 134 135 !------------ 136 ! diagnostics 137 !------------ 391 138 DO jj = 1, jpj 392 139 DO ji = 1, jpi … … 404 151 IF( iom_use('destrp') ) CALL iom_put( "destrp" , diag_trp_es ) ! advected snw enthalpy (W/m2) 405 152 153 !-------------------------------------- 154 ! Thickness correction in case too high 155 !-------------------------------------- 406 156 IF( nn_limdyn == 2 ) THEN 407 157 ! 408 CALL ice_var_zapsmall !--- zap small areas ---!158 CALL ice_var_zapsmall !-- zap small areas 409 159 ! 410 DO jl = 1, jpl !--- Thickness correction in case too high --- !160 DO jl = 1, jpl 411 161 DO jj = 1, jpj 412 162 DO ji = 1, jpi 413 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 163 IF ( v_i(ji,jj,jl) > 0._wp ) THEN !-- bound to zhimax 414 164 ! 415 165 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) … … 428 178 END DO 429 179 430 WHERE( ht_i(:,:,jpl) > hi_max(jpl) ) !--- bound ht_i to hi_max (99 m)180 WHERE( ht_i(:,:,jpl) > hi_max(jpl) ) !-- bound ht_i to hi_max (99 m) 431 181 ht_i(:,:,jpl) = hi_max(jpl) 432 182 a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl) 433 183 END WHERE 434 184 435 IF ( nn_pnd_scheme > 0 ) THEN !--- correct pond fraction to avoid a_ip > a_i185 IF ( nn_pnd_scheme > 0 ) THEN !-- correct pond fraction to avoid a_ip > a_i 436 186 WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:) 437 187 ENDIF … … 442 192 ! Impose a_i < amax if no ridging/rafting or in mono-category 443 193 !------------------------------------------------------------ 444 ! 445 IF( l_piling ) THEN ! simple conservative piling, comparable with 1-cat models 194 IF( l_piling ) THEN !-- simple conservative piling, comparable with 1-cat models 446 195 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 447 196 DO jl = 1, jpl … … 452 201 ENDIF 453 202 454 ! --- agglomerate variables -----------------203 ! agglomerate variables 455 204 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 456 205 vt_s(:,:) = SUM( v_s(:,:,:), dim=3 ) … … 464 213 ! END MP 2016 465 214 466 ! --- open water = 1 if at_i=0 --------------------------------215 ! open water = 1 if at_i=0 467 216 WHERE( at_i == 0._wp ) ato_i = 1._wp 468 217 … … 470 219 IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 471 220 472 ! -------------- -----------------------------------221 ! -------------- 473 222 ! control prints 474 ! -------------- -----------------------------------223 ! -------------- 475 224 IF( ln_limctl ) CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 476 225 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90
r8486 r8504 12 12 !! 'key_lim3' LIM3 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 !! ice_adv_x : advection of sea ice on x axis 15 !! ice_adv_y : advection of sea ice on y axis 14 !! ice_adv_prather : advection of sea ice using Prather scheme 16 15 !!---------------------------------------------------------------------- 17 16 USE dom_oce ! ocean domain 18 17 USE ice ! sea-ice variables 18 USE sbc_oce , ONLY : nn_fsbc ! frequency of sea-ice call 19 19 ! 20 20 USE lbclnk ! lateral boundary condition - MPP exchanges … … 27 27 PRIVATE 28 28 29 PUBLIC ice_adv_x ! called by iceadv 30 PUBLIC ice_adv_y ! called by iceadv 29 PUBLIC ice_adv_prather ! called by iceadv 31 30 32 31 !! * Substitutions … … 39 38 CONTAINS 40 39 41 SUBROUTINE ice_adv_ x( pdf, put , pcrh, psm , ps0 ,&42 & psx, psxx, psy , psyy, psxy)40 SUBROUTINE ice_adv_prather( kt, pu_ice, pv_ice, & 41 & pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 43 42 !!---------------------------------------------------------------------- 44 !! ** routine ice_adv_ x**43 !! ** routine ice_adv_prather ** 45 44 !! 46 45 !! ** purpose : Computes and adds the advection trend to sea-ice 47 !! variable on i-axis48 46 !! 49 47 !! ** method : Uses Prather second order scheme that advects tracers 50 !! but also theirquadratic forms. The method preserves51 !! tracer structures by conserving second order moments.48 !! but also their quadratic forms. The method preserves 49 !! tracer structures by conserving second order moments. 52 50 !! 53 51 !! Reference: Prather, 1986, JGR, 91, D6. 6671-6681. 54 52 !!---------------------------------------------------------------------- 53 INTEGER , INTENT(in ) :: kt ! time step 54 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity 55 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity 56 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 57 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume 58 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume 59 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psmv_i ! salt content 60 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content 61 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration 62 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 63 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 64 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 65 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 66 ! 67 INTEGER :: jk, jl, jt ! dummy loop indices 68 INTEGER :: initad ! number of sub-timestep for the advection 69 REAL(wp) :: zcfl , zusnit ! - - 70 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zarea 71 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0opw 72 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0es , z0smi , z0oi 73 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ap , z0vp 74 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z0ei 75 !!---------------------------------------------------------------------- 76 ! 77 IF( kt == nit000 .AND. lwp ) THEN 78 WRITE(numout,*) 79 WRITE(numout,*) 'ice_adv_prather: Prather advection scheme' 80 WRITE(numout,*) '~~~~~~~~~~~~~~~' 81 ENDIF 82 ! 83 ALLOCATE( zarea(jpi,jpj) , z0opw(jpi,jpj, 1 ) , & 84 & z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) , & 85 & z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) , & 86 & z0ei (jpi,jpj,nlay_i,jpl) ) 87 ! 88 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 89 zcfl = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 90 zcfl = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 91 IF( lk_mpp ) CALL mpp_max( zcfl ) 92 93 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 94 ELSE ; initad = 1 ; zusnit = 1.0_wp 95 ENDIF 96 97 zarea(:,:) = e1e2t(:,:) 98 !------------------------- 99 ! transported fields 100 !------------------------- 101 z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:) ! Open water area 102 DO jl = 1, jpl 103 z0snw(:,:,jl) = pv_s (:,:, jl) * e1e2t(:,:) ! Snow volume 104 z0ice(:,:,jl) = pv_i (:,:, jl) * e1e2t(:,:) ! Ice volume 105 z0ai (:,:,jl) = pa_i (:,:, jl) * e1e2t(:,:) ! Ice area 106 z0smi(:,:,jl) = psmv_i(:,:, jl) * e1e2t(:,:) ! Salt content 107 z0oi (:,:,jl) = poa_i (:,:, jl) * e1e2t(:,:) ! Age content 108 z0es (:,:,jl) = pe_s (:,:,1,jl) * e1e2t(:,:) ! Snow heat content 109 DO jk = 1, nlay_i 110 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 111 END DO 112 IF ( nn_pnd_scheme > 0 ) THEN 113 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 114 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 115 ENDIF 116 END DO 117 118 ! !--------------------------------------------! 119 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 120 ! !--------------------------------------------! 121 DO jt = 1, initad 122 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 123 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 124 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), & 125 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 126 DO jl = 1, jpl 127 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 128 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 129 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & 130 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 131 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 132 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 133 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & 134 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 135 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 136 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 137 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & 138 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 139 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 140 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 141 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & 142 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 143 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 144 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 145 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & 146 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 147 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 148 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 149 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & 150 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 151 DO jk = 1, nlay_i !--- ice heat contents --- 152 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 153 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 154 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 155 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 156 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 157 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 158 END DO 159 IF ( nn_pnd_scheme > 0 ) THEN 160 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & !--- melt pond fraction -- 161 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 162 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & 163 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 164 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & !--- melt pond volume -- 165 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 166 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & 167 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 168 ENDIF 169 END DO 170 END DO 171 ! !--------------------------------------------! 172 ELSE !== even ice time step: adv_y then adv_x ==! 173 ! !--------------------------------------------! 174 DO jt = 1, initad 175 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 176 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 177 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), & 178 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 179 DO jl = 1, jpl 180 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 181 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 182 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & 183 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 184 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 185 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 186 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & 187 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 188 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 189 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 190 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & 191 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 192 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 193 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 194 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & 195 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 196 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 197 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 198 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & 199 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 200 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 201 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 202 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & 203 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 204 DO jk = 1, nlay_i !--- ice heat contents --- 205 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 206 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 207 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 208 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 209 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 210 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 211 END DO 212 IF ( nn_pnd_scheme > 0 ) THEN 213 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & !--- melt pond fraction --- 214 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 215 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & 216 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) ) 217 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & !--- melt pond volume --- 218 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 219 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & 220 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) ) 221 ENDIF 222 END DO 223 END DO 224 ENDIF 225 226 !------------------------------------------- 227 ! Recover the properties from their contents 228 !------------------------------------------- 229 pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) 230 DO jl = 1, jpl 231 pv_i (:,:, jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) 232 pv_s (:,:, jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) 233 psmv_i(:,:, jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) 234 poa_i (:,:, jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) 235 pa_i (:,:, jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) 236 pe_s (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:) 237 DO jk = 1, nlay_i 238 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) 239 END DO 240 ! MV MP 2016 241 IF ( nn_pnd_scheme > 0 ) THEN 242 pa_ip (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) 243 pv_ip (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) 244 ENDIF 245 ! END MV MP 2016 246 END DO 247 ! 248 DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei ) 249 ! 250 END SUBROUTINE ice_adv_prather 251 252 SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 , & 253 & psx, psxx, psy , psyy, psxy ) 254 !!---------------------------------------------------------------------- 255 !! ** routine adv_x ** 256 !! 257 !! ** purpose : Computes and adds the advection trend to sea-ice 258 !! variable on x axis 259 !!---------------------------------------------------------------------- 55 260 REAL(wp) , INTENT(in ) :: pdf ! reduction factor for the time step 56 REAL(wp) , INTENT(in ) :: pcrh ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0)261 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 57 262 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s] 58 263 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psm ! area … … 209 414 210 415 IF(ln_ctl) THEN 211 CALL prt_ctl(tab2d_1=psm , clinfo1=' ice_adv_x: psm :', tab2d_2=ps0 , clinfo2=' ps0 : ')212 CALL prt_ctl(tab2d_1=psx , clinfo1=' ice_adv_x: psx :', tab2d_2=psxx, clinfo2=' psxx : ')213 CALL prt_ctl(tab2d_1=psy , clinfo1=' ice_adv_x: psy :', tab2d_2=psyy, clinfo2=' psyy : ')214 CALL prt_ctl(tab2d_1=psxy , clinfo1=' ice_adv_x: psxy :')416 CALL prt_ctl(tab2d_1=psm , clinfo1=' adv_x: psm :', tab2d_2=ps0 , clinfo2=' ps0 : ') 417 CALL prt_ctl(tab2d_1=psx , clinfo1=' adv_x: psx :', tab2d_2=psxx, clinfo2=' psxx : ') 418 CALL prt_ctl(tab2d_1=psy , clinfo1=' adv_x: psy :', tab2d_2=psyy, clinfo2=' psyy : ') 419 CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :') 215 420 ENDIF 216 421 ! 217 END SUBROUTINE ice_adv_x218 219 220 SUBROUTINE ice_adv_y( pdf, pvt , pcrh, psm , ps0 , &221 & 422 END SUBROUTINE adv_x 423 424 425 SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 , & 426 & psx, psxx, psy , psyy, psxy ) 222 427 !!--------------------------------------------------------------------- 223 !! ** routine ice_adv_y **428 !! ** routine adv_y ** 224 429 !! 225 430 !! ** purpose : Computes and adds the advection trend to sea-ice 226 !! variable on y axis 227 !! 228 !! ** method : Uses Prather second order scheme that advects tracers 229 !! but also their quadratic forms. The method preserves 230 !! tracer structures by conserving second order moments. 231 !! 232 !! Reference: Prather, 1986, JGR, 91, D6. 6671-6681. 431 !! variable on y axis 233 432 !!--------------------------------------------------------------------- 234 433 REAL(wp) , INTENT(in ) :: pdf ! reduction factor for the time step 235 REAL(wp) , INTENT(in ) :: pcrh ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0)434 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 236 435 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s] 237 436 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psm ! area … … 389 588 390 589 IF(ln_ctl) THEN 391 CALL prt_ctl(tab2d_1=psm , clinfo1=' ice_adv_y: psm :', tab2d_2=ps0 , clinfo2=' ps0 : ')392 CALL prt_ctl(tab2d_1=psx , clinfo1=' ice_adv_y: psx :', tab2d_2=psxx, clinfo2=' psxx : ')393 CALL prt_ctl(tab2d_1=psy , clinfo1=' ice_adv_y: psy :', tab2d_2=psyy, clinfo2=' psyy : ')394 CALL prt_ctl(tab2d_1=psxy , clinfo1=' ice_adv_y: psxy :')590 CALL prt_ctl(tab2d_1=psm , clinfo1=' adv_y: psm :', tab2d_2=ps0 , clinfo2=' ps0 : ') 591 CALL prt_ctl(tab2d_1=psx , clinfo1=' adv_y: psx :', tab2d_2=psxx, clinfo2=' psxx : ') 592 CALL prt_ctl(tab2d_1=psy , clinfo1=' adv_y: psy :', tab2d_2=psyy, clinfo2=' psyy : ') 593 CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :') 395 594 ENDIF 396 595 ! 397 END SUBROUTINE ice_adv_y596 END SUBROUTINE adv_y 398 597 399 598 #else -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90
r8486 r8504 43 43 CONTAINS 44 44 45 SUBROUTINE ice_adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc ) 45 SUBROUTINE ice_adv_umx( kt, pu_ice, pv_ice, & 46 & pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 46 47 !!---------------------------------------------------------------------- 47 48 !! *** ROUTINE ice_adv_umx *** 49 !! 50 !! ** Purpose : Compute the now trend due to total advection of 51 !! tracers and add it to the general trend of tracer equations 52 !! using an "Ultimate-Macho" scheme 53 !! 54 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 55 !!---------------------------------------------------------------------- 56 INTEGER , INTENT(in ) :: kt ! time step 57 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity 58 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity 59 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 60 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume 61 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume 62 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psmv_i ! salt content 63 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content 64 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration 65 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 66 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 67 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 68 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 69 ! 70 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 71 INTEGER :: initad ! number of sub-timestep for the advection 72 REAL(wp) :: zcfl , zusnit, zdt ! - - 73 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zudy, zvdx, zcu_box, zcv_box 74 !!---------------------------------------------------------------------- 75 ! 76 IF( kt == nit000 .AND. lwp ) THEN 77 WRITE(numout,*) 78 WRITE(numout,*) 'ice_adv_umx : Ultimate-MACHO advection scheme' 79 WRITE(numout,*) '~~~~~~~~~~~' 80 ENDIF 81 ! 82 ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 83 ! 84 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 85 zcfl = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 86 zcfl = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 87 IF( lk_mpp ) CALL mpp_max( zcfl ) 88 89 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 90 ELSE ; initad = 1 ; zusnit = 1.0_wp 91 ENDIF 92 93 zdt = rdt_ice / REAL(initad) 94 95 ! --- transport --- ! 96 zudy(:,:) = pu_ice(:,:) * e2u(:,:) 97 zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 98 99 ! --- define velocity for advection: u*grad(H) --- ! 100 DO jj = 2, jpjm1 101 DO ji = fs_2, fs_jpim1 102 IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp 103 ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj) 104 ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj) 105 ENDIF 106 107 IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp 108 ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1) 109 ELSE ; zcv_box(ji,jj) = pv_ice(ji,jj ) 110 ENDIF 111 END DO 112 END DO 113 114 !---------------! 115 !== advection ==! 116 !---------------! 117 DO jt = 1, initad 118 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) ) ! Open water area 119 DO jl = 1, jpl 120 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) ) ! Ice area 121 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) ) ! Ice volume 122 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) ) ! Salt content 123 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) ) ! Age content 124 DO jk = 1, nlay_i 125 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) ) ! Ice heat content 126 END DO 127 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) ) ! Snow volume 128 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) ) ! Snow heat content 129 IF ( nn_pnd_scheme > 0 ) THEN 130 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) ) ! Melt pond fraction 131 CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) ) ! Melt pond volume 132 ENDIF 133 END DO 134 END DO 135 ! 136 DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 137 ! 138 END SUBROUTINE ice_adv_umx 139 140 SUBROUTINE adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc ) 141 !!---------------------------------------------------------------------- 142 !! *** ROUTINE adv_umx *** 48 143 !! 49 144 !! ** Purpose : Compute the now trend due to total advection of … … 146 241 IF( nn_timing == 1 ) CALL timing_stop('ice_adv_umx') 147 242 ! 148 END SUBROUTINE ice_adv_umx243 END SUBROUTINE adv_umx 149 244 150 245
Note: See TracChangeset
for help on using the changeset viewer.
