Changeset 7646 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r6490 r7646 17 17 USE dom_oce ! ocean domain 18 18 USE sbc_oce ! ocean surface boundary condition 19 USE dom_ice ! ice domain20 19 USE ice ! ice variables 21 USE limadv ! ice advection22 20 USE limhdf ! ice horizontal diffusion 23 21 USE limvar ! 22 USE limadv_prather ! advection scheme (Prather) 23 USE limadv_umx ! advection scheme (ultimate-macho) 24 24 ! 25 25 USE in_out_manager ! I/O manager … … 57 57 !! ** method : variables included in the process are scalar, 58 58 !! other values are considered as second order. 59 !! For advection, a second order Prather scheme is used. 59 !! For advection, one can choose between 60 !! a) an Ultimate-Macho scheme (whose order is defined by nn_limadv_ord) => nn_limadv=0 61 !! b) and a second order Prather scheme => nn_limadv=-1 60 62 !! 61 63 !! ** action : 62 64 !!--------------------------------------------------------------------- 63 INTEGER, INTENT(in) :: kt 65 INTEGER, INTENT(in) :: kt ! number of iteration 64 66 ! 65 INTEGER :: ji, jj, jk, jm , jl, jt! dummy loop indices67 INTEGER :: ji, jj, jk, jm, jl, jt ! dummy loop indices 66 68 INTEGER :: initad ! number of sub-timestep for the advection 67 69 REAL(wp) :: zcfl , zusnit ! - - 68 CHARACTER(len=80) :: 70 CHARACTER(len=80) :: cltmp 69 71 ! 70 REAL(wp), POINTER, DIMENSION(:,:) :: zsm 72 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 73 REAL(wp) :: zdv, zda 74 REAL(wp), POINTER, DIMENSION(:,:) :: zatold, zeiold, zesold, zsmvold 75 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhimax, zviold, zvsold 76 ! --- diffusion --- ! 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhdfptab 78 INTEGER , PARAMETER :: ihdf_vars = 6 ! Number of variables in which we apply horizontal diffusion 79 ! inside limtrp for each ice category , not counting the 80 ! variables corresponding to ice_layers 81 ! --- ultimate macho only --- ! 82 REAL(wp) :: zdt 83 REAL(wp), POINTER, DIMENSION(:,:) :: zudy, zvdx, zcu_box, zcv_box 84 ! --- prather only --- ! 85 REAL(wp), POINTER, DIMENSION(:,:) :: zarea 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: z0opw 71 87 REAL(wp), POINTER, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0es , z0smi , z0oi 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: z0opw73 88 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: z0ei 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 75 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhimax ! old ice thickness 76 REAL(wp), POINTER, DIMENSION(:,:) :: zatold, zeiold, zesold ! old concentration, enthalpies 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhdfptab 78 REAL(wp) :: zdv, zvi, zvs, zsmv, zes, zei 79 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 80 !!--------------------------------------------------------------------- 81 INTEGER :: ihdf_vars = 6 !!Number of variables in which we apply horizontal diffusion 82 !! inside limtrp for each ice category , not counting the 83 !! variables corresponding to ice_layers 89 !! 84 90 !!--------------------------------------------------------------------- 85 91 IF( nn_timing == 1 ) CALL timing_start('limtrp') 86 92 87 CALL wrk_alloc( jpi,jpj, zsm, zatold, zeiold, zesold ) 88 CALL wrk_alloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 89 CALL wrk_alloc( jpi,jpj,1, z0opw ) 90 CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei ) 91 CALL wrk_alloc( jpi,jpj,jpl, zhimax, zviold, zvsold, zsmvold ) 92 CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1,zhdfptab) 93 94 IF( numit == nstart .AND. lwp ) THEN 95 WRITE(numout,*) 96 IF( ln_limdyn ) THEN ; WRITE(numout,*) 'lim_trp : Ice transport ' 97 ELSE ; WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn 98 ENDIF 99 WRITE(numout,*) '~~~~~~~~~~~~' 93 CALL wrk_alloc( jpi,jpj, zatold, zeiold, zesold, zsmvold ) 94 CALL wrk_alloc( jpi,jpj,jpl, zhimax, zviold, zvsold ) 95 CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1, zhdfptab) 96 97 IF( kt == nit000 .AND. lwp ) THEN 98 WRITE(numout,*)'' 99 WRITE(numout,*)'limtrp' 100 WRITE(numout,*)'~~~~~~' 100 101 ncfl = 0 ! nb of time step with CFL > 1/2 101 102 ENDIF 102 103 zsm(:,:) = e1e2t(:,:) 104 105 ! !-------------------------------------! 106 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! 107 ! !-------------------------------------! 108 109 ! conservation test 110 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 111 112 ! mass and salt flux init 113 zviold(:,:,:) = v_i(:,:,:) 114 zvsold(:,:,:) = v_s(:,:,:) 115 zsmvold(:,:,:) = smv_i(:,:,:) 116 zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 117 zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 118 119 !--- Thickness correction init. ------------------------------- 120 zatold(:,:) = SUM( a_i(:,:,:), dim=3 ) 121 DO jl = 1, jpl 122 DO jj = 1, jpj 123 DO ji = 1, jpi 124 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 125 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 126 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 103 104 CALL lim_var_agg( 1 ) ! integrated values + ato_i 105 106 !-------------------------------------! 107 ! Advection of sea ice properties ! 108 !-------------------------------------! 109 110 ! conservation test 111 IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 112 113 ! store old values for diag 114 zviold = v_i 115 zvsold = v_s 116 zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 ) 117 zeiold (:,:) = et_i 118 zesold (:,:) = et_s 119 120 !--- Thickness correction init. --- ! 121 zatold(:,:) = at_i 122 DO jl = 1, jpl 123 DO jj = 1, jpj 124 DO ji = 1, jpi 125 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 126 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 127 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 128 END DO 129 END DO 130 END DO 131 ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- ! 132 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 133 DO jl = 1, jpl 134 DO jj = 2, jpjm1 135 DO ji = 2, jpim1 136 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) 137 END DO 138 END DO 139 CALL lbc_lnk(zhimax(:,:,jl),'T',1.) 140 END DO 141 142 ! --- If ice drift field is too fast, use an appropriate time step for advection --- ! 143 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability 144 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 145 IF( lk_mpp ) CALL mpp_max( zcfl ) 146 147 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 148 ELSE ; initad = 1 ; zusnit = 1.0_wp 149 ENDIF 150 151 !! IF( zcfl > 0.5_wp .AND. lwp ) THEN 152 !! ncfl = ncfl + 1 153 !! IF( ncfl > 0 ) THEN 154 !! WRITE(cltmp,'(i6.1)') ncfl 155 !! CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 156 !! ENDIF 157 !! ENDIF 158 159 SELECT CASE ( nn_limadv ) 160 161 !=============================! 162 CASE ( 0 ) !== Ultimate-MACHO scheme ==! 163 !=============================! 164 165 CALL wrk_alloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box ) 166 167 IF( kt == nit000 .AND. lwp ) THEN 168 WRITE(numout,*)'' 169 WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme' 170 WRITE(numout,*)'~~~~~~~~~~~' 171 ENDIF 172 ! 173 zdt = rdt_ice / REAL(initad) 174 175 ! transport 176 zudy(:,:) = u_ice(:,:) * e2u(:,:) 177 zvdx(:,:) = v_ice(:,:) * e1v(:,:) 178 179 ! define velocity for advection: u*grad(H) 180 DO jj = 2, jpjm1 181 DO ji = fs_2, fs_jpim1 182 IF ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp 183 ELSEIF( u_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = u_ice(ji-1,jj) 184 ELSE ; zcu_box(ji,jj) = u_ice(ji ,jj) 185 ENDIF 186 187 IF ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp 188 ELSEIF( v_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = v_ice(ji,jj-1) 189 ELSE ; zcv_box(ji,jj) = v_ice(ji,jj ) 190 ENDIF 191 END DO 192 END DO 193 194 ! advection 195 DO jt = 1, initad 196 CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) ) ! Open water area 197 DO jl = 1, jpl 198 CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) ) ! Ice area 199 CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) ) ! Ice volume 200 CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) ) ! Salt content 201 CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) ) ! Age content 202 DO jk = 1, nlay_i 203 CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) ) ! Ice heat content 127 204 END DO 128 END DO 129 END DO 130 !--------------------------------------------------------------------- 131 ! Record max of the surrounding ice thicknesses for correction 132 ! in case advection creates ice too thick. 133 !--------------------------------------------------------------------- 134 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 135 DO jl = 1, jpl 136 DO jj = 2, jpjm1 137 DO ji = 2, jpim1 138 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) 139 END DO 140 END DO 141 CALL lbc_lnk(zhimax(:,:,jl),'T',1.) 142 END DO 143 144 !=============================! 145 !== Prather scheme ==! 146 !=============================! 147 148 ! If ice drift field is too fast, use an appropriate time step for advection. 149 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability 150 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 151 IF(lk_mpp ) CALL mpp_max( zcfl ) 152 153 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 154 ELSE ; initad = 1 ; zusnit = 1.0_wp 205 CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) ) ! Snow volume 206 CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) ) ! Snow heat content 207 END DO 208 END DO 209 ! 210 at_i(:,:) = a_i(:,:,1) ! total ice fraction 211 DO jl = 2, jpl 212 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 213 END DO 214 ! 215 CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box ) 216 217 !=============================! 218 CASE ( -1 ) !== Prather scheme ==! 219 !=============================! 220 221 CALL wrk_alloc( jpi,jpj, zarea ) 222 CALL wrk_alloc( jpi,jpj,1, z0opw ) 223 CALL wrk_alloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 224 CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei ) 225 226 IF( kt == nit000 .AND. lwp ) THEN 227 WRITE(numout,*)'' 228 WRITE(numout,*)'lim_adv_xy : Prather advection scheme' 229 WRITE(numout,*)'~~~~~~~~~~~' 155 230 ENDIF 156 157 IF( zcfl > 0.5_wp .AND. lwp ) ncfl = ncfl + 1 158 !! IF( lwp ) THEN 159 !! IF( ncfl > 0 ) THEN 160 !! WRITE(cltmp,'(i6.1)') ncfl 161 !! CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 162 !! ELSE 163 !! ! WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 ' 164 !! ENDIF 165 !! ENDIF 166 231 232 zarea(:,:) = e1e2t(:,:) 233 167 234 !------------------------- 168 235 ! transported fields … … 176 243 z0oi (:,:,jl) = oa_i (:,:, jl) * e1e2t(:,:) ! Age content 177 244 z0es (:,:,jl) = e_s (:,:,1,jl) * e1e2t(:,:) ! Snow heat content 178 DO jk = 1, nlay_i245 DO jk = 1, nlay_i 179 246 z0ei (:,:,jk,jl) = e_i (:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 180 247 END DO … … 184 251 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 185 252 DO jt = 1, initad 186 CALL lim_adv_x( zusnit, u_ice, 1._wp, z sm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area187 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )188 CALL lim_adv_y( zusnit, v_ice, 0._wp, z sm, z0opw (:,:,1), sxopw(:,:), &189 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )253 CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 254 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 255 CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), & 256 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 190 257 DO jl = 1, jpl 191 CALL lim_adv_x( zusnit, u_ice, 1._wp, z sm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---192 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )193 CALL lim_adv_y( zusnit, v_ice, 0._wp, z sm, z0ice (:,:,jl), sxice(:,:,jl), &194 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )195 CALL lim_adv_x( zusnit, u_ice, 1._wp, z sm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---196 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )197 CALL lim_adv_y( zusnit, v_ice, 0._wp, z sm, z0snw (:,:,jl), sxsn (:,:,jl), &198 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )199 CALL lim_adv_x( zusnit, u_ice, 1._wp, z sm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---200 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )201 CALL lim_adv_y( zusnit, v_ice, 0._wp, z sm, z0smi (:,:,jl), sxsal(:,:,jl), &202 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )203 CALL lim_adv_x( zusnit, u_ice, 1._wp, z sm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---204 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )205 CALL lim_adv_y( zusnit, v_ice, 0._wp, z sm, z0oi (:,:,jl), sxage(:,:,jl), &206 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )207 CALL lim_adv_x( zusnit, u_ice, 1._wp, z sm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---208 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )209 CALL lim_adv_y( zusnit, v_ice, 0._wp, z sm, z0ai (:,:,jl), sxa (:,:,jl), &210 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )211 CALL lim_adv_x( zusnit, u_ice, 1._wp, z sm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---212 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )213 CALL lim_adv_y( zusnit, v_ice, 0._wp, z sm, z0es (:,:,jl), sxc0 (:,:,jl), &214 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )258 CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 259 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 260 CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & 261 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 262 CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 263 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 264 CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & 265 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 266 CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 267 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 268 CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & 269 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 270 CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 271 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 272 CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & 273 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 274 CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 275 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 276 CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & 277 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 278 CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 279 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 280 CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & 281 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 215 282 DO jk = 1, nlay_i !--- ice heat contents --- 216 CALL lim_adv_x( zusnit, u_ice, 1._wp, z sm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &217 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), &218 & syye(:,:,jk,jl), sxye(:,:,jk,jl) )219 CALL lim_adv_y( zusnit, v_ice, 0._wp, z sm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &220 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), &221 & syye(:,:,jk,jl), sxye(:,:,jk,jl) )283 CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 284 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 285 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 286 CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 287 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 288 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 222 289 END DO 223 290 END DO … … 225 292 ELSE 226 293 DO jt = 1, initad 227 CALL lim_adv_y( zusnit, v_ice, 1._wp, z sm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area228 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )229 CALL lim_adv_x( zusnit, u_ice, 0._wp, z sm, z0opw (:,:,1), sxopw(:,:), &230 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )294 CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area 295 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 296 CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), & 297 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 231 298 DO jl = 1, jpl 232 CALL lim_adv_y( zusnit, v_ice, 1._wp, z sm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---233 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )234 CALL lim_adv_x( zusnit, u_ice, 0._wp, z sm, z0ice (:,:,jl), sxice(:,:,jl), &235 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )236 CALL lim_adv_y( zusnit, v_ice, 1._wp, z sm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---237 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )238 CALL lim_adv_x( zusnit, u_ice, 0._wp, z sm, z0snw (:,:,jl), sxsn (:,:,jl), &239 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )240 CALL lim_adv_y( zusnit, v_ice, 1._wp, z sm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---241 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )242 CALL lim_adv_x( zusnit, u_ice, 0._wp, z sm, z0smi (:,:,jl), sxsal(:,:,jl), &243 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )244 CALL lim_adv_y( zusnit, v_ice, 1._wp, z sm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---245 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )246 CALL lim_adv_x( zusnit, u_ice, 0._wp, z sm, z0oi (:,:,jl), sxage(:,:,jl), &247 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )248 CALL lim_adv_y( zusnit, v_ice, 1._wp, z sm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---249 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )250 CALL lim_adv_x( zusnit, u_ice, 0._wp, z sm, z0ai (:,:,jl), sxa (:,:,jl), &251 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )252 CALL lim_adv_y( zusnit, v_ice, 1._wp, z sm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---253 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )254 CALL lim_adv_x( zusnit, u_ice, 0._wp, z sm, z0es (:,:,jl), sxc0 (:,:,jl), &255 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )299 CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume --- 300 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 301 CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & 302 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 303 CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 304 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 305 CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & 306 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 307 CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 308 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 309 CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & 310 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 311 CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 312 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 313 CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & 314 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 315 CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 316 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 317 CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & 318 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 319 CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 320 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 321 CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es (:,:,jl), sxc0 (:,:,jl), & 322 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 256 323 DO jk = 1, nlay_i !--- ice heat contents --- 257 CALL lim_adv_y( zusnit, v_ice, 1._wp, z sm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &258 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), &259 & syye(:,:,jk,jl), sxye(:,:,jk,jl) )260 CALL lim_adv_x( zusnit, u_ice, 0._wp, z sm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &261 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), &262 & syye(:,:,jk,jl), sxye(:,:,jk,jl) )324 CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 325 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 326 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 327 CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), & 328 & sxxe(:,:,jk,jl), sye (:,:,jk,jl), & 329 & syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 263 330 END DO 264 331 END DO … … 286 353 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 287 354 END DO 288 289 !------------------------------------------------------------------------------! 290 ! Diffusion of Ice fields 291 !------------------------------------------------------------------------------! 292 !------------------------------------ 293 ! Diffusion of other ice variables 294 !------------------------------------ 355 356 CALL wrk_dealloc( jpi,jpj, zarea ) 357 CALL wrk_dealloc( jpi,jpj,1, z0opw ) 358 CALL wrk_dealloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 359 CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei ) 360 361 END SELECT 362 363 !------------------------------! 364 ! Diffusion of Ice fields 365 !------------------------------! 366 IF( nn_ahi0 /= -1 .AND. nn_limdyn == 2 ) THEN 367 ! 368 ! --- Prepare diffusion for variables with categories --- ! 369 ! mask eddy diffusivity coefficient at ocean U- and V-points 295 370 jm=1 296 371 DO jl = 1, jpl 297 ! ! Masked eddy diffusivity coefficient at ocean U- and V-points298 ! DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row299 ! DO ji = 1 , fs_jpim1 ! vector opt.300 ! pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji ,jj,jl) ) ) ) &301 ! & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)302 ! pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj ,jl) ) ) ) &303 ! & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)304 ! END DO305 ! END DO306 372 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 307 DO ji = 1 , fs_jpim1 ! vector opt.373 DO ji = 1 , fs_jpim1 308 374 pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji ,jj, jl ) ) ) ) & 309 375 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj, jl ) ) ) ) * ahiu(ji,jj) … … 313 379 END DO 314 380 315 zhdfptab(:,:,jm)= a_i (:,:, jl); jm = jm + 1 381 zhdfptab(:,:,jm)= a_i (:,:, jl); jm = jm + 1 316 382 zhdfptab(:,:,jm)= v_i (:,:, jl); jm = jm + 1 317 zhdfptab(:,:,jm)= v_s (:,:, jl); jm = jm + 1 383 zhdfptab(:,:,jm)= v_s (:,:, jl); jm = jm + 1 318 384 zhdfptab(:,:,jm)= smv_i(:,:, jl); jm = jm + 1 319 385 zhdfptab(:,:,jm)= oa_i (:,:, jl); jm = jm + 1 320 386 zhdfptab(:,:,jm)= e_s (:,:,1,jl); jm = jm + 1 321 ! Sample of adding more variables to apply lim_hdf using lim_hdf optimization--- 322 ! zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1 323 ! zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1 324 ! 325 ! and in this example the parameter ihdf_vars musb be changed to 8 (necessary for allocation) 326 !---------------------------------------------------------------------------------------- 387 ! Sample of adding more variables to apply lim_hdf (ihdf_vars must be increased) 388 ! zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1 389 ! zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1 327 390 DO jk = 1, nlay_i 328 391 zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1 329 392 END DO 330 393 END DO 331 ! 332 !-------------------------------- 333 ! diffusion of open water area 334 !-------------------------------- 335 ! ! Masked eddy diffusivity coefficient at ocean U- and V-points 336 !DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 337 ! DO ji = 1 , fs_jpim1 ! vector opt. 338 ! pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji ,jj) ) ) ) & 339 ! & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) 340 ! pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj ) ) ) ) & 341 ! & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj) 342 ! END DO 343 !END DO 344 394 395 ! --- Prepare diffusion for open water area --- ! 396 ! mask eddy diffusivity coefficient at ocean U- and V-points 345 397 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 346 DO ji = 1 , fs_jpim1 ! vector opt.398 DO ji = 1 , fs_jpim1 347 399 pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji ,jj) ) ) ) & 348 400 & * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) … … 353 405 ! 354 406 zhdfptab(:,:,jm)= ato_i (:,:); 355 CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i) 356 407 408 ! --- Apply diffusion --- ! 409 CALL lim_hdf( zhdfptab, ihdf_vars ) 410 411 ! --- Recover properties --- ! 357 412 jm=1 358 413 DO jl = 1, jpl 359 a_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 360 v_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 361 v_s (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 362 smv_i(:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 363 oa_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 364 e_s (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1 365 ! Sample of adding more variables to apply lim_hdf--------- 366 ! variable_1 (:,:,1,jl) = zhdfptab(:,:, jm ) ; jm + 1 367 ! variable_2 (:,:,1,jl) = zhdfptab(:,:, jm ) ; jm + 1 368 !----------------------------------------------------------- 414 a_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 415 v_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 416 v_s (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 417 smv_i(:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 418 oa_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1 419 e_s (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1 420 ! Sample of adding more variables to apply lim_hdf 421 ! variable_1 (:,:,1,jl) = zhdfptab(:,:, jm ) ; jm + 1 422 ! variable_2 (:,:,1,jl) = zhdfptab(:,:, jm ) ; jm + 1 369 423 DO jk = 1, nlay_i 370 e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 371 END DO 372 END DO 373 424 e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 425 END DO 426 END DO 374 427 ato_i (:,:) = zhdfptab(:,:,jm) 375 376 !------------------------------------------------------------------------------! 377 ! limit ice properties after transport 378 !------------------------------------------------------------------------------! 379 !!gm & cr : MAX should not be active if adv scheme is positive ! 428 429 ENDIF 430 431 ! --- diags --- 432 DO jj = 1, jpj 433 DO ji = 1, jpi 434 diag_trp_ei (ji,jj) = ( SUM( e_i (ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 435 diag_trp_es (ji,jj) = ( SUM( e_s (ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 436 diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:) ) - zsmvold(ji,jj) ) * r1_rdtice 437 diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 438 diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 439 END DO 440 END DO 441 442 IF( nn_limdyn == 2) THEN 443 444 ! zap small areas 445 CALL lim_var_zapsmall 446 447 !--- Thickness correction in case too high --- ! 380 448 DO jl = 1, jpl 381 449 DO jj = 1, jpj 382 450 DO ji = 1, jpi 383 v_s (ji,jj,jl) = MAX( 0._wp, v_s (ji,jj,jl) ) 384 v_i (ji,jj,jl) = MAX( 0._wp, v_i (ji,jj,jl) ) 385 smv_i(ji,jj,jl) = MAX( 0._wp, smv_i(ji,jj,jl) ) 386 oa_i (ji,jj,jl) = MAX( 0._wp, oa_i (ji,jj,jl) ) 387 a_i (ji,jj,jl) = MAX( 0._wp, a_i (ji,jj,jl) ) 388 e_s (ji,jj,1,jl) = MAX( 0._wp, e_s (ji,jj,1,jl) ) 389 END DO 390 END DO 391 392 DO jk = 1, nlay_i 393 DO jj = 1, jpj 394 DO ji = 1, jpi 395 e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) ) 396 END DO 397 END DO 398 END DO 399 END DO 400 !!gm & cr 401 402 ! --- diags --- 403 DO jj = 1, jpj 404 DO ji = 1, jpi 405 diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 406 diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 407 408 diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 409 diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 410 diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice 411 END DO 412 END DO 413 414 ! zap small areas 415 CALL lim_var_zapsmall 416 417 !--- Thickness correction in case too high -------------------------------------------------------- 418 DO jl = 1, jpl 419 DO jj = 1, jpj 420 DO ji = 1, jpi 421 451 422 452 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 423 453 424 454 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 425 455 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 426 456 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 427 457 428 zvi = v_i (ji,jj,jl)429 zvs = v_s (ji,jj,jl)430 zsmv = smv_i(ji,jj,jl)431 zes = e_s (ji,jj,1,jl)432 zei = SUM( e_i(ji,jj,1:nlay_i,jl) )433 434 458 zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 435 459 436 460 IF ( ( zdv > 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 437 461 & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 438 462 439 463 rswitch = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) 440 464 a_i(ji,jj,jl) = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 ) 441 465 442 466 ! small correction due to *rswitch for a_i 443 467 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl) … … 446 470 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 447 471 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 448 449 ! Update mass fluxes 450 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 451 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 452 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 453 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 454 hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 455 472 456 473 ENDIF 457 474 458 475 ENDIF 459 476 460 477 END DO 461 478 END DO … … 463 480 ! ------------------------------------------------- 464 481 465 !-------------------------------------- 466 ! Impose a_i < amax in mono-category 467 !-------------------------------------- 468 ! 469 IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2 470 DO jj = 1, jpj 471 DO ji = 1, jpi 472 a_i(ji,jj,1) = MIN( a_i(ji,jj,1), rn_amax_2d(ji,jj) ) 473 END DO 474 END DO 475 ENDIF 476 477 ! --- agglomerate variables ----------------- 478 vt_i (:,:) = 0._wp 479 vt_s (:,:) = 0._wp 480 at_i (:,:) = 0._wp 482 ! Force the upper limit of ht_i to always be < hi_max (99 m). 483 DO jj = 1, jpj 484 DO ji = 1, jpi 485 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) ) 486 ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) ) 487 a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch 488 END DO 489 END DO 490 491 ENDIF 492 493 !------------------------------------------------------------ 494 ! Impose a_i < amax if no ridging/rafting or in mono-category 495 !------------------------------------------------------------ 496 ! 497 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 498 IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2 481 499 DO jl = 1, jpl 482 500 DO jj = 1, jpj 483 501 DO ji = 1, jpi 484 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) 485 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) 486 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) 502 rswitch = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) ) 503 zda = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp ) & 504 & * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 ) 505 a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda 487 506 END DO 488 507 END DO 489 508 END DO 490 491 ! --- open water = 1 if at_i=0 --------------------------------492 DO jj = 1, jpj493 DO ji = 1, jpi494 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )495 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)496 END DO497 END DO498 499 ! conservation test500 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)501 502 509 ENDIF 503 510 511 ! --- agglomerate variables ----------------- 512 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 513 vt_s(:,:) = SUM( v_s(:,:,:), dim=3 ) 514 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 515 516 ! --- open water = 1 if at_i=0 -------------------------------- 517 WHERE( at_i == 0._wp ) ato_i = 1._wp 518 519 ! conservation test 520 IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 521 504 522 ! ------------------------------------------------- 505 523 ! control prints 506 524 ! ------------------------------------------------- 507 IF( ln_ icectl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )525 IF( ln_limctl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 508 526 ! 509 CALL wrk_dealloc( jpi,jpj, zsm, zatold, zeiold, zesold ) 510 CALL wrk_dealloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 511 CALL wrk_dealloc( jpi,jpj,1, z0opw ) 512 CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei ) 513 CALL wrk_dealloc( jpi,jpj,jpl, zviold, zvsold, zhimax, zsmvold ) 514 CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars+nlay_i)+1,zhdfptab) 527 CALL wrk_dealloc( jpi,jpj, zatold, zeiold, zesold, zsmvold ) 528 CALL wrk_dealloc( jpi,jpj,jpl, zhimax, zviold, zvsold ) 529 CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1, zhdfptab) 515 530 ! 516 531 IF( nn_timing == 1 ) CALL timing_stop('limtrp') 517 532 ! 518 533 END SUBROUTINE lim_trp 519 534
Note: See TracChangeset
for help on using the changeset viewer.
