Changeset 888 for trunk/NEMO/LIM_SRC_2/limrhg_2.F90
- Timestamp:
- 2008-04-11T19:05:03+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_2/limrhg_2.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_2/limrhg_2.F90
r823 r888 4 4 !! Ice rheology : performs sea ice rheology 5 5 !!====================================================================== 6 !! History : 0.0 ! 93-12 (M.A. Morales Maqueda.) Original code 7 !! 1.0 ! 94-12 (H. Goosse) 8 !! 2.0 ! 03-12 (C. Ethe, G. Madec) F90, mpp 9 !! " " ! 06-08 (G. Madec) surface module, ice-stress at I-point 10 !! " " ! 09-09 (G. Madec) Huge verctor optimisation 11 !!---------------------------------------------------------------------- 6 12 #if defined key_lim2 7 13 !!---------------------------------------------------------------------- 8 14 !! 'key_lim2' LIM 2.0 sea-ice model 9 15 !!---------------------------------------------------------------------- 16 !!---------------------------------------------------------------------- 10 17 !! lim_rhg_2 : computes ice velocities 11 18 !!---------------------------------------------------------------------- 12 !! * Modules used13 USE phycst14 USE par_oce15 USE ice_oce !ice variables16 USE dom_ice_217 USE ice_2 18 USE lbclnk 19 USE lib_mpp 20 USE in_out_manager 21 USE prtctl 19 USE par_oce ! ocean parameter 20 USE ice_oce ! ice variables 21 USE sbc_ice ! surface boundary condition: ice variables 22 USE dom_ice_2 ! domaine: ice variables 23 USE phycst ! physical constant 24 USE ice_2 ! ice variables 25 USE lbclnk ! lateral boundary condition 26 USE lib_mpp ! MPP library 27 USE in_out_manager ! I/O manager 28 USE prtctl ! Print control 22 29 23 30 IMPLICIT NONE 24 31 PRIVATE 25 32 26 !! * Routine accessibility27 PUBLIC lim_rhg_2 ! routine called by lim_dyn_2 28 29 !! * Module variables30 REAL(wp) :: & ! constant values 31 rzero = 0.e0 , &32 rone = 1.e0 33 !!---------------------------------------------------------------------- 34 !! LIM 2.0, UCL-LOCEAN-IPSL (200 5)35 !! $ Header$36 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt33 PUBLIC lim_rhg_2 ! routine called by lim_dyn 34 35 REAL(wp) :: rzero = 0.e0 ! constant value: zero 36 REAL(wp) :: rone = 1.e0 ! and one 37 38 !! * Substitutions 39 # include "vectopt_loop_substitute.h90" 40 !!---------------------------------------------------------------------- 41 !! LIM 2.0, UCL-LOCEAN-IPSL (2006) 42 !! $ Id: $ 43 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 37 44 !!---------------------------------------------------------------------- 38 45 … … 48 55 !! viscous-plastic law including shear strength and a bulk rheology. 49 56 !! 50 !! ** Action : - compute u_ice, v_ice the sea-ice velocity 57 !! ** Action : - compute ui_ice, vi_ice the sea-ice velocity defined 58 !! at I-point 59 !!------------------------------------------------------------------- 60 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation 61 INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation 51 62 !! 52 !! History : 53 !! 0.0 ! 93-12 (M.A. Morales Maqueda.) Original code 54 !! 1.0 ! 94-12 (H. Goosse) 55 !! 2.0 ! 03-12 (C. Ethe, G. Madec) F90, mpp 63 INTEGER :: ji, jj ! dummy loop indices 64 INTEGER :: iter, jter ! temporary integers 65 CHARACTER (len=50) :: charout 66 REAL(wp) :: ze11 , ze12 , ze22 , ze21 ! temporary scalars 67 REAL(wp) :: zt11 , zt12 , zt21 , zt22 ! " " 68 REAL(wp) :: zvis11, zvis21, zvis12, zvis22 ! " " 69 REAL(wp) :: zgphsx, ztagnx, zunw, zur, zusw ! " " 70 REAL(wp) :: zgphsy, ztagny, zvnw, zvr ! " " 71 REAL(wp) :: zresm, za, zac, zmod 72 REAL(wp) :: zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal 73 REAL(wp) :: ztrace2, zdeter, zdelta, zmask, zdgp, zdgi, zdiag 74 REAL(wp) :: za1, zb1, zc1, zd1 75 REAL(wp) :: za2, zb2, zc2, zd2, zden 76 REAL(wp) :: zs11_11, zs11_12, zs11_21, zs11_22 77 REAL(wp) :: zs12_11, zs12_12, zs12_21, zs12_22 78 REAL(wp) :: zs21_11, zs21_12, zs21_21, zs21_22 79 REAL(wp) :: zs22_11, zs22_12, zs22_21, zs22_22 80 REAL(wp), DIMENSION(jpi, jpj ) :: zfrld, zmass, zcorl 81 REAL(wp), DIMENSION(jpi, jpj ) :: za1ct, za2ct, zresr 82 REAL(wp), DIMENSION(jpi, jpj ) :: zc1u, zc1v, zc2u, zc2v 83 REAL(wp), DIMENSION(jpi, jpj ) :: zsang 84 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zu0, zv0 85 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zu_n, zv_n 86 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zu_a, zv_a 87 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zviszeta, zviseta 88 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zzfrld, zztms 89 REAL(wp), DIMENSION(jpi,0:jpj+1) :: zi1, zi2, zmasst, zpresh 90 56 91 !!------------------------------------------------------------------- 57 ! * Arguments 58 INTEGER, INTENT(in) :: k_j1 , & ! southern j-index for ice computation 59 & k_jpj ! northern j-index for ice computation 60 61 ! * Local variables 62 INTEGER :: ji, jj ! dummy loop indices 63 64 INTEGER :: & 65 iim1, ijm1, iip1 , ijp1 , & ! temporary integers 66 iter, jter ! " " 67 68 CHARACTER (len=50) :: charout 69 70 REAL(wp) :: & 71 ze11 , ze12 , ze22 , ze21 , & ! temporary scalars 72 zt11 , zt12 , zt21 , zt22 , & ! " " 73 zvis11, zvis21, zvis12, zvis22, & ! " " 74 zgphsx, ztagnx, zusw , & ! " " 75 zgphsy, ztagny ! " " 76 REAL(wp) :: & 77 zresm, zunw, zvnw, zur, zvr, zmod, za, zac, & 78 zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal, & 79 ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag 80 REAL(wp),DIMENSION(jpi,jpj) :: & 81 zpresh, zfrld, zmass, zcorl, & 82 zu0, zv0, zviszeta, zviseta, & 83 zc1u, zc1v, zc2u, zc2v, za1ct, za2ct, za1, za2, zb1, zb2, & 84 zc1, zc2, zd1, zd2, zden, zu_ice, zv_ice, zresr 85 REAL(wp),DIMENSION(jpi,jpj,2,2) :: & 86 zs11, zs12, zs22, zs21 87 !!------------------------------------------------------------------- 92 93 !!bug 94 !! ui_oce(:,:) = 0.e0 95 !! vi_oce(:,:) = 0.e0 96 !! write(*,*) 'rhg min, max u & v', maxval(ui_oce), minval(ui_oce), maxval(vi_oce), minval(vi_oce) 97 !!bug 88 98 89 99 ! Store initial velocities 90 ! ------------------------ 91 zu0(:,:) = u_ice(:,:) 92 zv0(:,:) = v_ice(:,:) 100 ! ---------------- 101 zztms(:,0 ) = 0.e0 ; zzfrld(:,0 ) = 0.e0 102 zztms(:,jpj+1) = 0.e0 ; zzfrld(:,jpj+1) = 0.e0 103 zu0(:,0 ) = 0.e0 ; zv0(:,0 ) = 0.e0 104 zu0(:,jpj+1) = 0.e0 ; zv0(:,jpj+1) = 0.e0 105 zztms(:,1:jpj) = tms(:,:) ; zzfrld(:,1:jpj) = frld(:,:) 106 zu0(:,1:jpj) = ui_ice(:,:) ; zv0(:,1:jpj) = vi_ice(:,:) 107 108 zu_a(:,:) = zu0(:,:) ; zv_a(:,:) = zv0(:,:) 109 zu_n(:,:) = zu0(:,:) ; zv_n(:,:) = zv0(:,:) 110 111 !i 112 zi1 (:,:) = 0.e0 113 zi2 (:,:) = 0.e0 114 zpresh(:,:) = 0.e0 115 zmasst(:,:) = 0.e0 116 !i 117 !!gm violant 118 zfrld(:,:) =0.e0 119 zcorl(:,:) =0.e0 120 zmass(:,:) =0.e0 121 za1ct(:,:) =0.e0 122 za2ct(:,:) =0.e0 123 !!gm end 124 125 zviszeta(:,:) = 0.e0 126 zviseta (:,:) = 0.e0 127 128 !i zviszeta(:,0 ) = 0.e0 ; zviseta(:,0 ) = 0.e0 129 !i zviszeta(:,jpj ) = 0.e0 ; zviseta(:,jpj ) = 0.e0 130 !i zviszeta(:,jpj+1) = 0.e0 ; zviseta(:,jpj+1) = 0.e0 131 93 132 94 133 ! Ice mass, ice strength, and wind stress at the center | … … 96 135 !------------------------------------------------------------------- 97 136 137 !CDIR NOVERRCHK 98 138 DO jj = k_j1 , k_jpj-1 139 !CDIR NOVERRCHK 99 140 DO ji = 1 , jpi 100 za1(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) 141 ! only the sinus changes its sign with the hemisphere 142 zsang(ji,jj) = SIGN( 1.e0, fcor(ji,jj) ) * sangvg ! only the sinus changes its sign with the hemisphere 143 ! 144 zmasst(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) 101 145 zpresh(ji,jj) = tms(ji,jj) * pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) ) 102 #if defined key_lim_cp1 && defined key_coupled 103 zb1(ji,jj) = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) ) 104 zb2(ji,jj) = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) ) 105 #else 106 zb1(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 107 zb2(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 108 #endif 146 !!gm :: stress given at I-point (F-point for the ocean) only compute the ponderation with the ice fraction (1-frld) 147 zi1(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 148 zi2(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 109 149 END DO 110 150 END DO … … 117 157 118 158 DO jj = k_j1+1, k_jpj-1 119 DO ji = 2, jpi120 zstms = tms(ji,jj ) * wght(ji,jj,2,2) +tms(ji-1,jj ) * wght(ji,jj,1,2) &121 & + tms(ji,jj-1) * wght(ji,jj,2,1) +tms(ji-1,jj-1) * wght(ji,jj,1,1)159 DO ji = fs_2, jpi 160 zstms = zztms(ji,jj ) * wght(ji,jj,2,2) + zztms(ji-1,jj ) * wght(ji,jj,1,2) & 161 & + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1) 122 162 zusw = 1.0 / MAX( zstms, epsd ) 123 163 124 zt11 = tms(ji ,jj ) *frld(ji ,jj )125 zt12 = tms(ji-1,jj ) *frld(ji-1,jj )126 zt21 = tms(ji ,jj-1) *frld(ji ,jj-1)127 zt22 = tms(ji-1,jj-1) *frld(ji-1,jj-1)164 zt11 = zztms(ji ,jj ) * zzfrld(ji ,jj ) 165 zt12 = zztms(ji-1,jj ) * zzfrld(ji-1,jj ) 166 zt21 = zztms(ji ,jj-1) * zzfrld(ji ,jj-1) 167 zt22 = zztms(ji-1,jj-1) * zzfrld(ji-1,jj-1) 128 168 129 169 ! Leads area. … … 131 171 & + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw 132 172 133 ! Mass and coriolis coeff. 134 zmass(ji,jj) = ( z a1(ji,jj ) * wght(ji,jj,2,2) + za1(ji-1,jj ) * wght(ji,jj,1,2) &135 & + z a1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw173 ! Mass and coriolis coeff. at I-point 174 zmass(ji,jj) = ( zmasst(ji,jj ) * wght(ji,jj,2,2) + zmasst(ji-1,jj ) * wght(ji,jj,1,2) & 175 & + zmasst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 136 176 zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj) 137 177 138 178 ! Wind stress. 139 #if defined key_lim_cp1 && defined key_coupled 140 ztagnx = ( zb1(ji,jj ) * wght(ji,jj,2,2) + zb1(ji-1,jj ) * wght(ji,jj,1,2) & 141 & + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 142 ztagny = ( zb2(ji,jj ) * wght(ji,jj,2,2) + zb2(ji-1,jj ) * wght(ji,jj,1,2) & 143 & + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 144 #else 145 ztagnx = ( zb1(ji,jj ) * wght(ji,jj,2,2) + zb1(ji-1,jj ) * wght(ji,jj,1,2) & 146 & + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj) 147 ztagny = ( zb2(ji,jj ) * wght(ji,jj,2,2) + zb2(ji-1,jj ) * wght(ji,jj,1,2) & 148 & + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj) 149 #endif 179 ! always provide stress at I-point (ocean F-point) 180 ztagnx = ( zi1(ji,jj ) * wght(ji,jj,2,2) + zi1(ji-1,jj ) * wght(ji,jj,1,2) & 181 & + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utaui_ice(ji,jj) 182 ztagny = ( zi2(ji,jj ) * wght(ji,jj,2,2) + zi2(ji-1,jj ) * wght(ji,jj,1,2) & 183 & + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtaui_ice(ji,jj) 150 184 151 185 ! Gradient of ice strength … … 161 195 162 196 ! Computation of the velocity field taking into account the ice-ice interaction. 163 ! Terms that are independent of the velocity field.164 za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v _oce(ji,jj) - zgphsx165 za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u _oce(ji,jj) - zgphsy197 ! Terms that are independent of the ice velocity field. 198 za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * vi_oce(ji,jj) - zgphsx 199 za2ct(ji,jj) = ztagny + zcorl(ji,jj) * ui_oce(ji,jj) - zgphsy 166 200 END DO 167 201 END DO 168 169 !! inutile!!170 !!?? CALL lbc_lnk( za1ct, 'I', -1. )171 !!?? CALL lbc_lnk( za2ct, 'I', -1. )172 202 173 203 … … 182 212 ! Computation of free drift field for free slip boundary conditions. 183 213 184 DO jj = k_j1, k_jpj-1 185 DO ji = 1, jpim1 186 !- Rate of strain tensor. 187 zt11 = akappa(ji,jj,1,1) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) - u_ice(ji,jj ) - u_ice(ji ,jj+1) ) & 188 & + akappa(ji,jj,1,2) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) + v_ice(ji,jj ) + v_ice(ji ,jj+1) ) 189 zt12 = - akappa(ji,jj,2,2) * ( u_ice(ji ,jj) + u_ice(ji+1,jj ) - u_ice(ji,jj+1) - u_ice(ji+1,jj+1) ) & 190 & - akappa(ji,jj,2,1) * ( v_ice(ji ,jj) + v_ice(ji+1,jj ) + v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) 191 zt22 = - akappa(ji,jj,2,2) * ( v_ice(ji ,jj) + v_ice(ji+1,jj ) - v_ice(ji,jj+1) - v_ice(ji+1,jj+1) ) & 192 & + akappa(ji,jj,2,1) * ( u_ice(ji ,jj) + u_ice(ji+1,jj ) + u_ice(ji,jj+1) + u_ice(ji+1,jj+1) ) 193 zt21 = akappa(ji,jj,1,1) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) - v_ice(ji,jj ) - v_ice(ji ,jj+1) ) & 194 & - akappa(ji,jj,1,2) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) + u_ice(ji,jj ) + u_ice(ji ,jj+1) ) 195 196 !- Rate of strain tensor. 197 zdgp = zt11 + zt22 198 zdgi = zt12 + zt21 199 ztrace2 = zdgp * zdgp 200 zdeter = zt11 * zt22 - 0.25 * zdgi * zdgi 201 202 ! Creep limit depends on the size of the grid. 203 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2), creepl) 204 205 !- Computation of viscosities. 206 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 207 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 208 END DO 209 END DO 210 !!?? CALL lbc_lnk( zviszeta, 'I', -1. ) ! or T point??? semble reellement inutile 211 !!?? CALL lbc_lnk( zviseta , 'I', -1. ) 212 213 214 !- Determination of zc1u, zc2u, zc1v and zc2v. 215 DO jj = k_j1+1, k_jpj-1 216 DO ji = 2, jpim1 217 ze11 = akappa(ji-1,jj-1,1,1) 218 ze12 = +akappa(ji-1,jj-1,2,2) 219 ze22 = akappa(ji-1,jj-1,2,1) 220 ze21 = -akappa(ji-1,jj-1,1,2) 221 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 222 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 223 zvis12 = zviseta (ji-1,jj-1) + dm 224 zvis21 = zviseta (ji-1,jj-1) 225 226 zdiag = zvis22 * ( ze11 + ze22 ) 227 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 228 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 229 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 230 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 231 232 ze11 = -akappa(ji,jj-1,1,1) 233 ze12 = +akappa(ji,jj-1,2,2) 234 ze22 = akappa(ji,jj-1,2,1) 235 ze21 = -akappa(ji,jj-1,1,2) 236 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 237 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 238 zvis12 = zviseta (ji,jj-1) + dm 239 zvis21 = zviseta (ji,jj-1) 240 241 zdiag = zvis22 * ( ze11 + ze22 ) 242 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 243 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 244 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 245 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 246 247 ze11 = akappa(ji-1,jj,1,1) 248 ze12 = -akappa(ji-1,jj,2,2) 249 ze22 = akappa(ji-1,jj,2,1) 250 ze21 = -akappa(ji-1,jj,1,2) 251 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 252 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 253 zvis12 = zviseta (ji-1,jj) + dm 254 zvis21 = zviseta (ji-1,jj) 255 256 zdiag = zvis22 * ( ze11 + ze22 ) 257 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 258 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 259 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 260 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 261 262 ze11 = -akappa(ji,jj,1,1) 263 ze12 = -akappa(ji,jj,2,2) 264 ze22 = akappa(ji,jj,2,1) 265 ze21 = -akappa(ji,jj,1,2) 266 zvis11 = 2.0 * zviseta (ji,jj) + dm 267 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 268 zvis12 = zviseta (ji,jj) + dm 269 zvis21 = zviseta (ji,jj) 270 271 zdiag = zvis22 * ( ze11 + ze22 ) 272 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 273 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 274 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 275 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 276 END DO 277 END DO 278 279 DO jj = k_j1+1, k_jpj-1 280 DO ji = 2, jpim1 281 zc1u(ji,jj) = & 282 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 283 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 284 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 285 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 286 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 287 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 288 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 289 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 290 291 zc2u(ji,jj) = & 292 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 293 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 294 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 295 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 296 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 297 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 298 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 299 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 300 END DO 301 END DO 302 303 DO jj = k_j1+1, k_jpj-1 304 DO ji = 2, jpim1 305 ! zc1v , zc2v. 306 ze11 = akappa(ji-1,jj-1,1,2) 307 ze12 = -akappa(ji-1,jj-1,2,1) 308 ze22 = +akappa(ji-1,jj-1,2,2) 309 ze21 = akappa(ji-1,jj-1,1,1) 310 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 311 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 312 zvis12 = zviseta (ji-1,jj-1) + dm 313 zvis21 = zviseta (ji-1,jj-1) 314 315 zdiag = zvis22 * ( ze11 + ze22 ) 316 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 317 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 318 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 319 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 214 !CDIR NOVERRCHK 215 DO jj = k_j1, k_jpj-1 216 !CDIR NOVERRCHK 217 DO ji = 1, fs_jpim1 218 !- Rate of strain tensor. 219 zt11 = akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji,jj ) - zu_a(ji ,jj+1) ) & 220 & + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji,jj ) + zv_a(ji ,jj+1) ) 221 zt12 = - akappa(ji,jj,2,2) * ( zu_a(ji ,jj) + zu_a(ji+1,jj ) - zu_a(ji,jj+1) - zu_a(ji+1,jj+1) ) & 222 & - akappa(ji,jj,2,1) * ( zv_a(ji ,jj) + zv_a(ji+1,jj ) + zv_a(ji,jj+1) + zv_a(ji+1,jj+1) ) 223 zt22 = - akappa(ji,jj,2,2) * ( zv_a(ji ,jj) + zv_a(ji+1,jj ) - zv_a(ji,jj+1) - zv_a(ji+1,jj+1) ) & 224 & + akappa(ji,jj,2,1) * ( zu_a(ji ,jj) + zu_a(ji+1,jj ) + zu_a(ji,jj+1) + zu_a(ji+1,jj+1) ) 225 zt21 = akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji,jj ) - zv_a(ji ,jj+1) ) & 226 & - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji,jj ) + zu_a(ji ,jj+1) ) 227 228 !- Rate of strain tensor. 229 zdgp = zt11 + zt22 230 zdgi = zt12 + zt21 231 ztrace2 = zdgp * zdgp 232 zdeter = zt11 * zt22 - 0.25 * zdgi * zdgi 233 234 ! Creep limit depends on the size of the grid. 235 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2 ), creepl) 236 237 !- Computation of viscosities. 238 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 239 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 240 END DO 241 END DO 242 243 !- Determination of zc1u, zc2u, zc1v and zc2v. 244 DO jj = k_j1+1, k_jpj-1 245 DO ji = fs_2, fs_jpim1 246 !* zc1u , zc2v 247 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 248 zvis12 = zviseta (ji-1,jj-1) + dm 249 zvis21 = zviseta (ji-1,jj-1) 250 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 251 zdiag = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) ) 252 zs11_11 = zvis11 * akappa(ji-1,jj-1,1,1) + zdiag 253 zs12_11 = zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2) 254 zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2) 255 zs22_11 = zvis11 * akappa(ji-1,jj-1,2,1) + zdiag 256 257 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 258 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 259 zvis12 = zviseta (ji,jj-1) + dm 260 zvis21 = zviseta (ji,jj-1) 261 zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) ) 262 zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag 263 zs12_21 = zvis12 * akappa(ji,jj-1,2,2) - zvis21 * akappa(ji,jj-1,1,2) 264 zs22_21 = zvis11 * akappa(ji,jj-1,2,1) + zdiag 265 zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2) 266 267 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 268 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 269 zvis12 = zviseta (ji-1,jj) + dm 270 zvis21 = zviseta (ji-1,jj) 271 zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) ) 272 zs11_12 = zvis11 * akappa(ji-1,jj,1,1) + zdiag 273 zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2) 274 zs22_12 = zvis11 * akappa(ji-1,jj,2,1) + zdiag 275 zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2) 276 277 zvis11 = 2.0 * zviseta (ji,jj) + dm 278 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 279 zvis12 = zviseta (ji,jj) + dm 280 zvis21 = zviseta (ji,jj) 281 zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) ) 282 zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag 283 zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2) 284 zs22_22 = zvis11 * akappa(ji,jj,2,1) + zdiag 285 zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2) 286 287 zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 288 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 289 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 290 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 291 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 292 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 293 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 294 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 295 296 zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 297 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 298 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 299 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 300 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 301 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 302 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 303 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 304 305 !* zc1v , zc2v. 306 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 307 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 308 zvis12 = zviseta (ji-1,jj-1) + dm 309 zvis21 = zviseta (ji-1,jj-1) 310 zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) ) 311 zs11_11 = zvis11 * akappa(ji-1,jj-1,1,2) + zdiag 312 zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1) 313 zs22_11 = zvis11 * akappa(ji-1,jj-1,2,2) + zdiag 314 zs21_11 = zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1) 320 315 321 ze11 = akappa(ji,jj-1,1,2) 322 ze12 = -akappa(ji,jj-1,2,1) 323 ze22 = +akappa(ji,jj-1,2,2) 324 ze21 = -akappa(ji,jj-1,1,1) 325 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 326 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 327 zvis12 = zviseta (ji,jj-1) + dm 328 zvis21 = zviseta (ji,jj-1) 329 330 zdiag = zvis22 * ( ze11 + ze22 ) 331 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 332 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 333 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 334 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 335 336 ze11 = akappa(ji-1,jj,1,2) 337 ze12 = -akappa(ji-1,jj,2,1) 338 ze22 = -akappa(ji-1,jj,2,2) 339 ze21 = akappa(ji-1,jj,1,1) 340 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 341 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 342 zvis12 = zviseta (ji-1,jj) + dm 343 zvis21 = zviseta (ji-1,jj) 344 345 zdiag = zvis22 * ( ze11 + ze22 ) 346 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 347 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 348 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 349 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 350 351 ze11 = akappa(ji,jj,1,2) 352 ze12 = -akappa(ji,jj,2,1) 353 ze22 = -akappa(ji,jj,2,2) 354 ze21 = -akappa(ji,jj,1,1) 355 zvis11 = 2.0 * zviseta (ji,jj) + dm 356 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 357 zvis12 = zviseta (ji,jj) + dm 358 zvis21 = zviseta (ji,jj) 359 360 zdiag = zvis22 * ( ze11 + ze22 ) 361 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 362 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 363 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 364 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 365 366 END DO 367 END DO 368 369 DO jj = k_j1+1, k_jpj-1 370 DO ji = 2, jpim1 371 zc1v(ji,jj) = & 372 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 373 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 374 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 375 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 376 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 377 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 378 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 379 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 380 zc2v(ji,jj) = & 381 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 382 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 383 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 384 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 385 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 386 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 387 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 388 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 389 END DO 390 END DO 391 392 ! Relaxation. 393 394 iflag: DO jter = 1 , nbitdr 395 396 ! Store previous drift field. 397 DO jj = k_j1, k_jpj-1 398 zu_ice(:,jj) = u_ice(:,jj) 399 zv_ice(:,jj) = v_ice(:,jj) 316 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 317 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 318 zvis12 = zviseta (ji,jj-1) + dm 319 zvis21 = zviseta (ji,jj-1) 320 zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) ) 321 zs11_21 = zvis11 * akappa(ji,jj-1,1,2) + zdiag 322 zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1) 323 zs22_21 = zvis11 * akappa(ji,jj-1,2,2) + zdiag 324 zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1) 325 326 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 327 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 328 zvis12 = zviseta (ji-1,jj) + dm 329 zvis21 = zviseta (ji-1,jj) 330 zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) ) 331 zs11_12 = zvis11 * akappa(ji-1,jj,1,2) + zdiag 332 zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1) 333 zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag 334 zs21_12 = zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1) 335 336 zvis11 = 2.0 * zviseta (ji,jj) + dm 337 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 338 zvis12 = zviseta (ji,jj) + dm 339 zvis21 = zviseta (ji,jj) 340 zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) ) 341 zs11_22 = zvis11 * akappa(ji,jj,1,2) + zdiag 342 zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1) 343 zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag 344 zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1) 345 346 zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 347 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 348 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 349 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 350 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 351 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 352 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 353 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 354 355 zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 356 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 357 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 358 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 359 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 360 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 361 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 362 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 400 363 END DO 401 364 END DO 365 366 ! GAUSS-SEIDEL method 367 ! ! ================ ! 368 iflag: DO jter = 1 , nbitdr ! Relaxation ! 369 ! ! ================ ! 370 !CDIR NOVERRCHK 402 371 DO jj = k_j1+1, k_jpj-1 403 zsang = SIGN( 1.e0, fcor(1,jj) ) * sangvg ! only the sinus changes its sign with the hemisphere 404 DO ji = 2, jpim1 405 zur = u_ice(ji,jj) - u_oce(ji,jj) 406 zvr = v_ice(ji,jj) - v_oce(ji,jj) 407 zmod = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) ) 408 za = rhoco * zmod 409 zac = za * cangvg 410 zmpzas = alpha * zcorl(ji,jj) + za * zsang 372 !CDIR NOVERRCHK 373 DO ji = fs_2, fs_jpim1 374 ! 375 ze11 = akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj) 376 ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj) 377 ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj) 378 ze21 = akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj) 379 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 380 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 381 zvis12 = zviseta (ji,jj-1) + dm 382 zvis21 = zviseta (ji,jj-1) 383 zdiag = zvis22 * ( ze11 + ze22 ) 384 zs11_21 = zvis11 * ze11 + zdiag 385 zs12_21 = zvis12 * ze12 + zvis21 * ze21 386 zs22_21 = zvis11 * ze22 + zdiag 387 zs21_21 = zvis12 * ze21 + zvis21 * ze12 388 389 ze11 = akappa(ji-1,jj,1,1) * ( zu_a(ji ,jj+1) - zu_a(ji-1,jj+1) ) & 390 & + akappa(ji-1,jj,1,2) * ( zv_a(ji ,jj+1) + zv_a(ji-1,jj+1) ) 391 ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) ) & 392 & - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) ) 393 ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) ) & 394 & + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) ) 395 ze21 = akappa(ji-1,jj,1,1) * ( zv_a(ji ,jj+1) - zv_a(ji-1,jj+1) ) & 396 & - akappa(ji-1,jj,1,2) * ( zu_a(ji ,jj+1) + zu_a(ji-1,jj+1) ) 397 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 398 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 399 zvis12 = zviseta (ji-1,jj) + dm 400 zvis21 = zviseta (ji-1,jj) 401 zdiag = zvis22 * ( ze11 + ze22 ) 402 zs11_12 = zvis11 * ze11 + zdiag 403 zs12_12 = zvis12 * ze12 + zvis21 * ze21 404 zs22_12 = zvis11 * ze22 + zdiag 405 zs21_12 = zvis12 * ze21 + zvis21 * ze12 406 407 ze11 = akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji ,jj+1) ) & 408 & + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji ,jj+1) ) 409 ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji ,jj+1) - zu_a(ji+1,jj+1) ) & 410 & - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji ,jj+1) + zv_a(ji+1,jj+1) ) 411 ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji ,jj+1) - zv_a(ji+1,jj+1) ) & 412 & + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji ,jj+1) + zu_a(ji+1,jj+1) ) 413 ze21 = akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji ,jj+1) ) & 414 & - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji ,jj+1) ) 415 zvis11 = 2.0 * zviseta (ji,jj) + dm 416 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 417 zvis12 = zviseta (ji,jj) + dm 418 zvis21 = zviseta (ji,jj) 419 zdiag = zvis22 * ( ze11 + ze22 ) 420 zs11_22 = zvis11 * ze11 + zdiag 421 zs12_22 = zvis12 * ze12 + zvis21 * ze21 422 zs22_22 = zvis11 * ze22 + zdiag 423 zs21_22 = zvis12 * ze21 + zvis21 * ze12 424 425 ! 2nd part 426 ze11 = akappa(ji-1,jj-1,1,1) * ( zu_a(ji ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) ) & 427 & + akappa(ji-1,jj-1,1,2) * ( zv_a(ji ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) ) 428 ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) - zu_a(ji-1,jj) ) & 429 & - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) + zv_a(ji-1,jj) ) 430 ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) - zv_a(ji-1,jj) ) & 431 & + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) + zu_a(ji-1,jj) ) 432 ze21 = akappa(ji-1,jj-1,1,1) * ( zv_a(ji ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) ) & 433 & - akappa(ji-1,jj-1,1,2) * ( zu_a(ji ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) ) 434 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 435 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 436 zvis12 = zviseta (ji-1,jj-1) + dm 437 zvis21 = zviseta (ji-1,jj-1) 438 zdiag = zvis22 * ( ze11 + ze22 ) 439 zs11_11 = zvis11 * ze11 + zdiag 440 zs12_11 = zvis12 * ze12 + zvis21 * ze21 441 zs22_11 = zvis11 * ze22 + zdiag 442 zs21_11 = zvis12 * ze21 + zvis21 * ze12 443 444 ze11 = akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji ,jj-1) ) & 445 & + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji ,jj-1) ) 446 ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) ) & 447 & - akappa(ji,jj-1,2,1) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) ) 448 ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) ) & 449 & + akappa(ji,jj-1,2,1) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) ) 450 ze21 = akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji ,jj-1) ) & 451 & - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji ,jj-1) ) 452 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 453 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 454 zvis12 = zviseta (ji,jj-1) + dm 455 zvis21 = zviseta (ji,jj-1) 456 zdiag = zvis22 * ( ze11 + ze22 ) 457 zs11_21 = zs11_21 + zvis11 * ze11 + zdiag 458 zs12_21 = zs12_21 + zvis12 * ze12 + zvis21 * ze21 459 zs22_21 = zs22_21 + zvis11 * ze22 + zdiag 460 zs21_21 = zs21_21 + zvis12 * ze21 + zvis21 * ze12 461 462 ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj) 463 ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj) 464 ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj) 465 ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj) 466 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 467 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 468 zvis12 = zviseta (ji-1,jj) + dm 469 zvis21 = zviseta (ji-1,jj) 470 zdiag = zvis22 * ( ze11 + ze22 ) 471 zs11_12 = zs11_12 + zvis11 * ze11 + zdiag 472 zs12_12 = zs12_12 + zvis12 * ze12 + zvis21 * ze21 473 zs22_12 = zs22_12 + zvis11 * ze22 + zdiag 474 zs21_12 = zs21_12 + zvis12 * ze21 + zvis21 * ze12 475 476 zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 477 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 478 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 479 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 480 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 481 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 482 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 483 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 484 485 zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 486 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 487 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 488 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 489 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 490 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 491 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 492 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 493 494 zur = zu_a(ji,jj) - ui_oce(ji,jj) 495 zvr = zv_a(ji,jj) - vi_oce(ji,jj) 496 !!!! 497 zmod = SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 498 za = rhoco * zmod 499 !!!! 500 !!gm chg resul za = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 501 zac = za * cangvg 502 zmpzas = alpha * zcorl(ji,jj) + za * zsang(ji,jj) 411 503 zmassdt = zusdtp * zmass(ji,jj) 412 504 zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj) 413 505 414 za1(ji,jj) = zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj) & 415 & + za * ( cangvg * u_oce(ji,jj) - zsang * v_oce(ji,jj) ) 416 417 za2(ji,jj) = zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj) & 418 & + za * ( cangvg * v_oce(ji,jj) + zsang * u_oce(ji,jj) ) 419 420 zb1(ji,jj) = zmassdt + zac - zc1u(ji,jj) 421 zb2(ji,jj) = zmpzas - zc2u(ji,jj) 422 zc1(ji,jj) = zmpzas + zc1v(ji,jj) 423 zc2(ji,jj) = zmassdt + zac - zc2v(ji,jj) 424 zdeter = zc1(ji,jj) * zb2(ji,jj) + zc2(ji,jj) * zb1(ji,jj) 425 zden(ji,jj) = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 506 za1 = zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj) & 507 & + za * ( cangvg * ui_oce(ji,jj) - zsang(ji,jj) * vi_oce(ji,jj) ) 508 za2 = zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj) & 509 & + za * ( cangvg * vi_oce(ji,jj) + zsang(ji,jj) * ui_oce(ji,jj) ) 510 zb1 = zmassdt + zac - zc1u(ji,jj) 511 zb2 = zmpzas - zc2u(ji,jj) 512 zc1 = zmpzas + zc1v(ji,jj) 513 zc2 = zmassdt + zac - zc2v(ji,jj) 514 zdeter = zc1 * zb2 + zc2 * zb1 515 zden = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 516 zunw = ( ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden 517 zvnw = ( ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden 518 zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 519 520 zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask 521 zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask 426 522 END DO 427 523 END DO 428 524 429 ! The computation of ice interaction term is splitted into two parts 430 !------------------------------------------------------------------------- 431 432 ! Terms that do not involve already up-dated velocities. 433 434 DO jj = k_j1+1, k_jpj-1 435 DO ji = 2, jpim1 436 iim1 = ji 437 ijm1 = jj - 1 438 iip1 = ji + 1 439 ijp1 = jj 440 ze11 = akappa(iim1,ijm1,1,1) * u_ice(iip1,ijp1) + akappa(iim1,ijm1,1,2) * v_ice(iip1,ijp1) 441 ze12 = + akappa(iim1,ijm1,2,2) * u_ice(iip1,ijp1) - akappa(iim1,ijm1,2,1) * v_ice(iip1,ijp1) 442 ze22 = + akappa(iim1,ijm1,2,2) * v_ice(iip1,ijp1) + akappa(iim1,ijm1,2,1) * u_ice(iip1,ijp1) 443 ze21 = akappa(iim1,ijm1,1,1) * v_ice(iip1,ijp1) - akappa(iim1,ijm1,1,2) * u_ice(iip1,ijp1) 444 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 445 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 446 zvis12 = zviseta (iim1,ijm1) + dm 447 zvis21 = zviseta (iim1,ijm1) 448 zdiag = zvis22 * ( ze11 + ze22 ) 449 zs11(ji,jj,2,1) = zvis11 * ze11 + zdiag 450 zs12(ji,jj,2,1) = zvis12 * ze12 + zvis21 * ze21 451 zs22(ji,jj,2,1) = zvis11 * ze22 + zdiag 452 zs21(ji,jj,2,1) = zvis12 * ze21 + zvis21 * ze12 453 454 455 iim1 = ji - 1 456 ijm1 = jj 457 iip1 = ji 458 ijp1 = jj + 1 459 ze11 = akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijp1) - u_ice(iim1,ijp1) ) & 460 & + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijp1) + v_ice(iim1,ijp1) ) 461 ze12 = + akappa(iim1,ijm1,2,2) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) & 462 & - akappa(iim1,ijm1,2,1) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) 463 ze22 = + akappa(iim1,ijm1,2,2) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) & 464 & + akappa(iim1,ijm1,2,1) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) 465 ze21 = akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijp1) - v_ice(iim1,ijp1) ) & 466 & - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 467 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 468 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 469 zvis12 = zviseta (iim1,ijm1) + dm 470 zvis21 = zviseta (iim1,ijm1) 471 zdiag = zvis22 * ( ze11 + ze22 ) 472 zs11(ji,jj,1,2) = zvis11 * ze11 + zdiag 473 zs12(ji,jj,1,2) = zvis12 * ze12 + zvis21 * ze21 474 zs22(ji,jj,1,2) = zvis11 * ze22 + zdiag 475 zs21(ji,jj,1,2) = zvis12 * ze21 + zvis21 * ze12 476 477 iim1 = ji 478 ijm1 = jj 479 iip1 = ji + 1 480 ijp1 = jj + 1 481 ze11 = akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) - u_ice(iim1,ijp1) ) & 482 & + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) + v_ice(iim1,ijp1) ) 483 ze12 = - akappa(iim1,ijm1,2,2) * ( u_ice(iip1,ijm1) - u_ice(iim1,ijp1) - u_ice(iip1,ijp1) ) & 484 & - akappa(iim1,ijm1,2,1) * ( v_ice(iip1,ijm1) + v_ice(iim1,ijp1) + v_ice(iip1,ijp1) ) 485 ze22 = - akappa(iim1,ijm1,2,2) * ( v_ice(iip1,ijm1) - v_ice(iim1,ijp1) - v_ice(iip1,ijp1) ) & 486 & + akappa(iim1,ijm1,2,1) * ( u_ice(iip1,ijm1) + u_ice(iim1,ijp1) + u_ice(iip1,ijp1) ) 487 ze21 = akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) - v_ice(iim1,ijp1) ) & 488 & - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 489 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 490 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 491 zvis12 = zviseta (iim1,ijm1) + dm 492 zvis21 = zviseta (iim1,ijm1) 493 494 zdiag = zvis22 * ( ze11 + ze22 ) 495 zs11(ji,jj,2,2) = zvis11 * ze11 + zdiag 496 zs12(ji,jj,2,2) = zvis12 * ze12 + zvis21 * ze21 497 zs22(ji,jj,2,2) = zvis11 * ze22 + zdiag 498 zs21(ji,jj,2,2) = zvis12 * ze21 + zvis21 * ze12 499 500 END DO 525 CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. ) 526 CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. ) 527 528 ! Test of Convergence 529 DO jj = k_j1+1 , k_jpj-1 530 zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) ) 501 531 END DO 502 503 ! Terms involving already up-dated velocities. 504 !-Using the arrays zu_ice and zv_ice in the computation of the terms ze leads to JACOBI's method; 505 ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method. 506 507 DO jj = k_j1+1, k_jpj-1 508 DO ji = 2, jpim1 509 iim1 = ji - 1 510 ijm1 = jj - 1 511 iip1 = ji 512 ijp1 = jj 513 ze11 = akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) - zu_ice(iim1,ijp1) ) & 514 & + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) + zv_ice(iim1,ijp1) ) 515 ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) - zu_ice(iim1,ijp1) ) & 516 & - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) + zv_ice(iim1,ijp1) ) 517 ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) - zv_ice(iim1,ijp1) ) & 518 & + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) + zu_ice(iim1,ijp1) ) 519 ze21 = akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) - zv_ice(iim1,ijp1) ) & 520 & - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) + zu_ice(iim1,ijp1) ) 521 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 522 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 523 zvis12 = zviseta (iim1,ijm1) + dm 524 zvis21 = zviseta (iim1,ijm1) 525 526 zdiag = zvis22 * ( ze11 + ze22 ) 527 zs11(ji,jj,1,1) = zvis11 * ze11 + zdiag 528 zs12(ji,jj,1,1) = zvis12 * ze12 + zvis21 * ze21 529 zs22(ji,jj,1,1) = zvis11 * ze22 + zdiag 530 zs21(ji,jj,1,1) = zvis12 * ze21 + zvis21 * ze12 531 532 #if defined key_agrif 533 END DO 534 END DO 535 536 DO jj = k_j1+1, k_jpj-1 537 DO ji = 2, jpim1 538 #endif 539 540 iim1 = ji 541 ijm1 = jj - 1 542 iip1 = ji + 1 543 ze11 = akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) ) & 544 & + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) ) 545 ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) ) & 546 & - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) ) 547 ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) ) & 548 & + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) ) 549 ze21 = akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) ) & 550 & - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) ) 551 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 552 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 553 zvis12 = zviseta (iim1,ijm1) + dm 554 zvis21 = zviseta (iim1,ijm1) 555 556 zdiag = zvis22 * ( ze11 + ze22 ) 557 zs11(ji,jj,2,1) = zs11(ji,jj,2,1) + zvis11 * ze11 + zdiag 558 zs12(ji,jj,2,1) = zs12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21 559 zs22(ji,jj,2,1) = zs22(ji,jj,2,1) + zvis11 * ze22 + zdiag 560 zs21(ji,jj,2,1) = zs21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12 561 562 563 iim1 = ji - 1 564 ijm1 = jj 565 ze11 = - akappa(iim1,ijm1,1,1) * zu_ice(iim1,ijm1) + akappa(iim1,ijm1,1,2) * zv_ice(iim1,ijm1) 566 ze12 = - akappa(iim1,ijm1,2,2) * zu_ice(iim1,ijm1) - akappa(iim1,ijm1,2,1) * zv_ice(iim1,ijm1) 567 ze22 = - akappa(iim1,ijm1,2,2) * zv_ice(iim1,ijm1) + akappa(iim1,ijm1,2,1) * zu_ice(iim1,ijm1) 568 ze21 = - akappa(iim1,ijm1,1,1) * zv_ice(iim1,ijm1) - akappa(iim1,ijm1,1,2) * zu_ice(iim1,ijm1) 569 zvis11 = 2.0 * zviseta (iim1,ijm1) + dm 570 zvis22 = zviszeta(iim1,ijm1) - zviseta(iim1,ijm1) 571 zvis12 = zviseta (iim1,ijm1) + dm 572 zvis21 = zviseta (iim1,ijm1) 573 574 zdiag = zvis22 * ( ze11 + ze22 ) 575 zs11(ji,jj,1,2) = zs11(ji,jj,1,2) + zvis11 * ze11 + zdiag 576 zs12(ji,jj,1,2) = zs12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21 577 zs22(ji,jj,1,2) = zs22(ji,jj,1,2) + zvis11 * ze22 + zdiag 578 zs21(ji,jj,1,2) = zs21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12 579 580 #if defined key_agrif 581 END DO 582 END DO 583 584 DO jj = k_j1+1, k_jpj-1 585 DO ji = 2, jpim1 586 #endif 587 zd1(ji,jj) = & 588 + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2) & 589 - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2) & 590 - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1) & 591 + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2) & 592 + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1) & 593 + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2) & 594 - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1) & 595 - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2) 596 zd2(ji,jj) = & 597 + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2) & 598 - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2) & 599 - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1) & 600 + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2) & 601 - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1) & 602 - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2) & 603 + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1) & 604 + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2) 605 END DO 532 zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) ) 533 !!!! this should be faster on scalar processor 534 ! zresm = MAXVAL( MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ), & 535 ! & ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) ) ) 536 !!!! 537 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 538 539 DO jj = k_j1, k_jpj 540 zu_a(:,jj) = zu_n(:,jj) 541 zv_a(:,jj) = zv_n(:,jj) 606 542 END DO 607 543 608 DO jj = k_j1+1, k_jpj-1 609 DO ji = 2, jpim1 610 zunw = ( ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj) & 611 & + ( za2(ji,jj) + zd2(ji,jj) ) * zc1(ji,jj) ) * zden(ji,jj) 612 613 zvnw = ( ( za2(ji,jj) + zd2(ji,jj) ) * zb1(ji,jj) & 614 & - ( za1(ji,jj) + zd1(ji,jj) ) * zb2(ji,jj) ) * zden(ji,jj) 615 616 zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 617 618 u_ice(ji,jj) = ( u_ice(ji,jj) + om * ( zunw - u_ice(ji,jj) ) * tmu(ji,jj) ) * zmask 619 v_ice(ji,jj) = ( v_ice(ji,jj) + om * ( zvnw - v_ice(ji,jj) ) * tmu(ji,jj) ) * zmask 620 END DO 544 IF( zresm <= resl ) EXIT iflag 545 546 ! ! ================ ! 547 END DO iflag ! end Relaxation ! 548 ! ! ================ ! 549 550 IF( zindu == 0 ) THEN ! even iteration 551 DO jj = k_j1 , k_jpj-1 552 zu0(:,jj) = zu_a(:,jj) 553 zv0(:,jj) = zv_a(:,jj) 621 554 END DO 622 623 CALL lbc_lnk( u_ice, 'I', -1. ) 624 CALL lbc_lnk( v_ice, 'I', -1. ) 625 626 !--- 5.2.5.4. Convergence test. 627 DO jj = k_j1+1 , k_jpj-1 628 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 629 END DO 630 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 631 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 632 633 IF ( zresm <= resl) EXIT iflag 634 635 END DO iflag 636 637 zindu1 = 1.0 - zindu 638 DO jj = k_j1 , k_jpj-1 639 zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj) 640 zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj) 641 END DO 642 ! ! ==================== ! 555 ENDIF 556 ! ! ==================== ! 643 557 END DO ! end loop over iter ! 644 558 ! ! ==================== ! 559 560 ui_ice(:,:) = zu_a(:,1:jpj) 561 vi_ice(:,:) = zv_a(:,1:jpj) 645 562 646 563 IF(ln_ctl) THEN 647 564 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter 648 565 CALL prt_ctl_info(charout) 649 CALL prt_ctl(tab2d_1=u _ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')566 CALL prt_ctl(tab2d_1=ui_ice, clinfo1=' lim_rhg : ui_ice :', tab2d_2=vi_ice, clinfo2=' vi_ice :') 650 567 ENDIF 651 568
Note: See TracChangeset
for help on using the changeset viewer.
