Changeset 719 for trunk/NEMO/LIM_SRC/limrhg.F90
- Timestamp:
- 2007-10-16T16:59:56+02:00 (17 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC/limrhg.F90 (modified) (7 diffs, 1 prop)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC/limrhg.F90
- Property svn:keywords changed from Id to Author Date Id Revision
r717 r719 4 4 !! Ice rheology : performs sea ice rheology 5 5 !!====================================================================== 6 !! History : 0.0 ! 93-12 (M.A. Morales Maqueda.) Original code7 !! 1.0 ! 94-12 (H. Goosse)8 !! 2.0 ! 03-12 (C. Ethe, G. Madec) F90, mpp9 !! " " ! 06-08 (G. Madec) surface module, ice-stress at I-point10 !! " " ! 09-09 (G. Madec) Huge verctor optimisation11 !!----------------------------------------------------------------------12 6 #if defined key_ice_lim 13 7 !!---------------------------------------------------------------------- 14 8 !! 'key_ice_lim' LIM sea-ice model 15 9 !!---------------------------------------------------------------------- 16 !!----------------------------------------------------------------------17 10 !! lim_rhg : computes ice velocities 18 11 !!---------------------------------------------------------------------- 19 USE par_oce ! ocean parameter20 USE ice_oce ! ice variables21 USE sbc_ice ! surface boundary condition: ice variables22 USE dom_ice ! domaine:ice variables23 USE phycst ! physical constant24 USE ice ! ice variables25 USE lbclnk ! lateral boundary condition26 USE lib_mpp ! MPP library27 USE in_out_manager ! I/O manager28 USE prtctl ! Print control12 !! * Modules used 13 USE phycst 14 USE par_oce 15 USE ice_oce ! ice variables 16 USE dom_ice 17 USE ice 18 USE lbclnk 19 USE lib_mpp 20 USE in_out_manager ! I/O manager 21 USE prtctl ! Print control 29 22 30 23 IMPLICIT NONE 31 24 PRIVATE 32 25 33 PUBLIC lim_rhg ! routine called by lim_dyn34 35 REAL(wp) :: rzero = 0.e0 ! constant value: zero 36 REAL(wp) :: rone = 1.e0 ! and one37 38 !! * Substitutions39 # include "vectopt_loop_substitute.h90" 26 !! * Routine accessibility 27 PUBLIC lim_rhg ! routine called by lim_dyn 28 29 !! * Module variables 30 REAL(wp) :: & ! constant values 31 rzero = 0.e0 , & 32 rone = 1.e0 40 33 !!---------------------------------------------------------------------- 41 !! LIM 2.0, UCL-LOCEAN-IPSL (200 6)42 !! $Header : $43 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)34 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 35 !! $Header$ 36 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 44 37 !!---------------------------------------------------------------------- 45 38 … … 55 48 !! viscous-plastic law including shear strength and a bulk rheology. 56 49 !! 57 !! ** Action : - compute ui_ice, vi_ice the sea-ice velocity defined 58 !! at I-point 50 !! ** Action : - compute u_ice, v_ice the sea-ice velocity 51 !! 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 59 56 !!------------------------------------------------------------------- 60 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation 61 INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation 62 !! 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 63 62 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 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 91 87 !!------------------------------------------------------------------- 92 93 !!bug94 !! ui_oce(:,:) = 0.e095 !! vi_oce(:,:) = 0.e096 !! write(*,*) 'rhg min, max u & v', maxval(ui_oce), minval(ui_oce), maxval(vi_oce), minval(vi_oce)97 !!bug98 88 99 89 ! Store initial velocities 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 90 ! ------------------------ 91 zu0(:,:) = u_ice(:,:) 92 zv0(:,:) = v_ice(:,:) 132 93 133 94 ! Ice mass, ice strength, and wind stress at the center | … … 135 96 !------------------------------------------------------------------- 136 97 137 !CDIR NOVERRCHK138 98 DO jj = k_j1 , k_jpj-1 139 !CDIR NOVERRCHK140 99 DO ji = 1 , jpi 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) ) 100 za1(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) ) 145 101 zpresh(ji,jj) = tms(ji,jj) * pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) ) 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) ) 149 150 !!gm #if defined key_lim_cp1 && defined key_coupled 151 !!gm zi1(ji,jj) = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) ) 152 !!gm zi2(ji,jj) = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) ) 153 !!gm #else 154 !!gm zi1(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 155 !!gm zi2(ji,jj) = tms(ji,jj) * ( 1.0 - frld(ji,jj) ) 156 !!gm #endif 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 157 109 END DO 158 110 END DO … … 165 117 166 118 DO jj = k_j1+1, k_jpj-1 167 DO ji = fs_2, jpi168 zstms = zztms(ji,jj ) * wght(ji,jj,2,2) + zztms(ji-1,jj ) * wght(ji,jj,1,2) &169 & + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1)119 DO ji = 2, jpi 120 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) 170 122 zusw = 1.0 / MAX( zstms, epsd ) 171 123 172 zt11 = zztms(ji ,jj ) * zzfrld(ji ,jj )173 zt12 = zztms(ji-1,jj ) * zzfrld(ji-1,jj )174 zt21 = zztms(ji ,jj-1) * zzfrld(ji ,jj-1)175 zt22 = zztms(ji-1,jj-1) * zzfrld(ji-1,jj-1)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) 176 128 177 129 ! Leads area. … … 179 131 & + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw 180 132 181 ! Mass and coriolis coeff. at I-point182 zmass(ji,jj) = ( z masst(ji,jj ) * wght(ji,jj,2,2) + zmasst(ji-1,jj ) * wght(ji,jj,1,2) &183 & + z masst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw133 ! Mass and coriolis coeff. 134 zmass(ji,jj) = ( za1(ji,jj ) * wght(ji,jj,2,2) + za1(ji-1,jj ) * wght(ji,jj,1,2) & 135 & + za1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 184 136 zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj) 185 137 186 138 ! Wind stress. 187 !!gm #if defined key_lim_cp1 && defined key_coupled 188 !!gm ztagnx = ( zi1(ji,jj ) * wght(ji,jj,2,2) + zi1(ji-1,jj ) * wght(ji,jj,1,2) & 189 !!gm & + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 190 !!gm ztagny = ( zi2(ji,jj ) * wght(ji,jj,2,2) + zi2(ji-1,jj ) * wght(ji,jj,1,2) & 191 !!gm & + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw 192 !!gm #else 193 !!gm ztagnx = ( zi1(ji,jj ) * wght(ji,jj,2,2) + zi1(ji-1,jj ) * wght(ji,jj,1,2) & 194 !!gm & + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj) 195 !!gm ztagny = ( zi2(ji,jj ) * wght(ji,jj,2,2) + zi2(ji-1,jj ) * wght(ji,jj,1,2) & 196 !!gm & + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj) 197 !!gm #endif 198 !!gm always stress provided at I-point (ocean F-point) 199 ztagnx = ( zi1(ji,jj ) * wght(ji,jj,2,2) + zi1(ji-1,jj ) * wght(ji,jj,1,2) & 200 & + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utaui_ice(ji,jj) 201 ztagny = ( zi2(ji,jj ) * wght(ji,jj,2,2) + zi2(ji-1,jj ) * wght(ji,jj,1,2) & 202 & + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtaui_ice(ji,jj) 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 203 150 204 151 ! Gradient of ice strength … … 214 161 215 162 ! Computation of the velocity field taking into account the ice-ice interaction. 216 ! Terms that are independent of the icevelocity field.217 za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v i_oce(ji,jj) - zgphsx218 za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u i_oce(ji,jj) - zgphsy163 ! Terms that are independent of the velocity field. 164 za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_oce(ji,jj) - zgphsx 165 za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_oce(ji,jj) - zgphsy 219 166 END DO 220 167 END DO 168 169 !! inutile!! 170 !!?? CALL lbc_lnk( za1ct, 'I', -1. ) 171 !!?? CALL lbc_lnk( za2ct, 'I', -1. ) 221 172 222 173 … … 231 182 ! Computation of free drift field for free slip boundary conditions. 232 183 233 !CDIR NOVERRCHK 234 DO jj = k_j1, k_jpj-1 235 !CDIR NOVERRCHK 236 DO ji = 1, fs_jpim1 237 !- Rate of strain tensor. 238 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) ) & 239 & + 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) ) 240 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) ) & 241 & - 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) ) 242 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) ) & 243 & + 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) ) 244 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) ) & 245 & - 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) ) 246 247 !- Rate of strain tensor. 248 zdgp = zt11 + zt22 249 zdgi = zt12 + zt21 250 ztrace2 = zdgp * zdgp 251 zdeter = zt11 * zt22 - 0.25 * zdgi * zdgi 252 253 ! Creep limit depends on the size of the grid. 254 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2 ), creepl) 255 256 !- Computation of viscosities. 257 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn ) 258 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2 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 320 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) 259 400 END DO 260 END DO 261 262 !- Determination of zc1u, zc2u, zc1v and zc2v. 263 DO jj = k_j1+1, k_jpj-1 264 DO ji = fs_2, fs_jpim1 265 !* zc1u , zc2v 266 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 267 zvis12 = zviseta (ji-1,jj-1) + dm 268 zvis21 = zviseta (ji-1,jj-1) 269 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 270 zdiag = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) ) 271 zs11_11 = zvis11 * akappa(ji-1,jj-1,1,1) + zdiag 272 zs12_11 = zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2) 273 zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2) 274 zs22_11 = zvis11 * akappa(ji-1,jj-1,2,1) + zdiag 275 276 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 277 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 278 zvis12 = zviseta (ji,jj-1) + dm 279 zvis21 = zviseta (ji,jj-1) 280 zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) ) 281 zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag 282 zs12_21 = zvis12 * akappa(ji,jj-1,2,2) - zvis21 * akappa(ji,jj-1,1,2) 283 zs22_21 = zvis11 * akappa(ji,jj-1,2,1) + zdiag 284 zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2) 285 286 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 287 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 288 zvis12 = zviseta (ji-1,jj) + dm 289 zvis21 = zviseta (ji-1,jj) 290 zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) ) 291 zs11_12 = zvis11 * akappa(ji-1,jj,1,1) + zdiag 292 zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2) 293 zs22_12 = zvis11 * akappa(ji-1,jj,2,1) + zdiag 294 zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2) 295 296 zvis11 = 2.0 * zviseta (ji,jj) + dm 297 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 298 zvis12 = zviseta (ji,jj) + dm 299 zvis21 = zviseta (ji,jj) 300 zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) ) 301 zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag 302 zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2) 303 zs22_22 = zvis11 * akappa(ji,jj,2,1) + zdiag 304 zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2) 305 306 zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 307 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 308 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 309 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 310 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 311 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 312 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 313 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 314 315 zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 316 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 317 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 318 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 319 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 320 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 321 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 322 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 323 324 !* zc1v , zc2v. 325 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 326 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 327 zvis12 = zviseta (ji-1,jj-1) + dm 328 zvis21 = zviseta (ji-1,jj-1) 329 zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) ) 330 zs11_11 = zvis11 * akappa(ji-1,jj-1,1,2) + zdiag 331 zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1) 332 zs22_11 = zvis11 * akappa(ji-1,jj-1,2,2) + zdiag 333 zs21_11 = zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1) 334 335 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 336 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 337 zvis12 = zviseta (ji,jj-1) + dm 338 zvis21 = zviseta (ji,jj-1) 339 zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) ) 340 zs11_21 = zvis11 * akappa(ji,jj-1,1,2) + zdiag 341 zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1) 342 zs22_21 = zvis11 * akappa(ji,jj-1,2,2) + zdiag 343 zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1) 344 345 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 346 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 347 zvis12 = zviseta (ji-1,jj) + dm 348 zvis21 = zviseta (ji-1,jj) 349 zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) ) 350 zs11_12 = zvis11 * akappa(ji-1,jj,1,2) + zdiag 351 zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1) 352 zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag 353 zs21_12 = zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1) 354 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 zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) ) 360 zs11_22 = zvis11 * akappa(ji,jj,1,2) + zdiag 361 zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1) 362 zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag 363 zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1) 364 365 zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 366 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 367 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 368 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 369 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 370 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 371 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 372 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 373 374 zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 375 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 376 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 377 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 378 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 379 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 380 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 381 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 382 END DO 383 END DO 384 385 ! GAUSS-SEIDEL method 386 ! ! ================ ! 387 iflag: DO jter = 1 , nbitdr ! Relaxation ! 388 ! ! ================ ! 389 !CDIR NOVERRCHK 401 390 402 DO jj = k_j1+1, k_jpj-1 391 !CDIR NOVERRCHK 392 DO ji = fs_2, fs_jpim1 393 ! 394 ze11 = akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj) 395 ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj) 396 ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj) 397 ze21 = akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj) 398 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 399 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 400 zvis12 = zviseta (ji,jj-1) + dm 401 zvis21 = zviseta (ji,jj-1) 402 zdiag = zvis22 * ( ze11 + ze22 ) 403 zs11_21 = zvis11 * ze11 + zdiag 404 zs12_21 = zvis12 * ze12 + zvis21 * ze21 405 zs22_21 = zvis11 * ze22 + zdiag 406 zs21_21 = zvis12 * ze21 + zvis21 * ze12 407 408 ze11 = akappa(ji-1,jj,1,1) * ( zu_a(ji ,jj+1) - zu_a(ji-1,jj+1) ) & 409 & + akappa(ji-1,jj,1,2) * ( zv_a(ji ,jj+1) + zv_a(ji-1,jj+1) ) 410 ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) ) & 411 & - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) ) 412 ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) ) & 413 & + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) ) 414 ze21 = akappa(ji-1,jj,1,1) * ( zv_a(ji ,jj+1) - zv_a(ji-1,jj+1) ) & 415 & - akappa(ji-1,jj,1,2) * ( zu_a(ji ,jj+1) + zu_a(ji-1,jj+1) ) 416 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 417 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 418 zvis12 = zviseta (ji-1,jj) + dm 419 zvis21 = zviseta (ji-1,jj) 420 zdiag = zvis22 * ( ze11 + ze22 ) 421 zs11_12 = zvis11 * ze11 + zdiag 422 zs12_12 = zvis12 * ze12 + zvis21 * ze21 423 zs22_12 = zvis11 * ze22 + zdiag 424 zs21_12 = zvis12 * ze21 + zvis21 * ze12 425 426 ze11 = akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji ,jj+1) ) & 427 & + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji ,jj+1) ) 428 ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji ,jj+1) - zu_a(ji+1,jj+1) ) & 429 & - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji ,jj+1) + zv_a(ji+1,jj+1) ) 430 ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji ,jj+1) - zv_a(ji+1,jj+1) ) & 431 & + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji ,jj+1) + zu_a(ji+1,jj+1) ) 432 ze21 = akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji ,jj+1) ) & 433 & - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji ,jj+1) ) 434 zvis11 = 2.0 * zviseta (ji,jj) + dm 435 zvis22 = zviszeta(ji,jj) - zviseta(ji,jj) 436 zvis12 = zviseta (ji,jj) + dm 437 zvis21 = zviseta (ji,jj) 438 zdiag = zvis22 * ( ze11 + ze22 ) 439 zs11_22 = zvis11 * ze11 + zdiag 440 zs12_22 = zvis12 * ze12 + zvis21 * ze21 441 zs22_22 = zvis11 * ze22 + zdiag 442 zs21_22 = zvis12 * ze21 + zvis21 * ze12 443 444 ! 2nd part 445 ze11 = akappa(ji-1,jj-1,1,1) * ( zu_a(ji ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) ) & 446 & + akappa(ji-1,jj-1,1,2) * ( zv_a(ji ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) ) 447 ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) - zu_a(ji-1,jj) ) & 448 & - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) + zv_a(ji-1,jj) ) 449 ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) - zv_a(ji-1,jj) ) & 450 & + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) + zu_a(ji-1,jj) ) 451 ze21 = akappa(ji-1,jj-1,1,1) * ( zv_a(ji ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) ) & 452 & - akappa(ji-1,jj-1,1,2) * ( zu_a(ji ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) ) 453 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm 454 zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1) 455 zvis12 = zviseta (ji-1,jj-1) + dm 456 zvis21 = zviseta (ji-1,jj-1) 457 zdiag = zvis22 * ( ze11 + ze22 ) 458 zs11_11 = zvis11 * ze11 + zdiag 459 zs12_11 = zvis12 * ze12 + zvis21 * ze21 460 zs22_11 = zvis11 * ze22 + zdiag 461 zs21_11 = zvis12 * ze21 + zvis21 * ze12 462 463 ze11 = akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji ,jj-1) ) & 464 & + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji ,jj-1) ) 465 ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) ) & 466 & - akappa(ji,jj-1,2,1) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) ) 467 ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) ) & 468 & + akappa(ji,jj-1,2,1) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) ) 469 ze21 = akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji ,jj-1) ) & 470 & - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji ,jj-1) ) 471 zvis11 = 2.0 * zviseta (ji,jj-1) + dm 472 zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1) 473 zvis12 = zviseta (ji,jj-1) + dm 474 zvis21 = zviseta (ji,jj-1) 475 zdiag = zvis22 * ( ze11 + ze22 ) 476 zs11_21 = zs11_21 + zvis11 * ze11 + zdiag 477 zs12_21 = zs12_21 + zvis12 * ze12 + zvis21 * ze21 478 zs22_21 = zs22_21 + zvis11 * ze22 + zdiag 479 zs21_21 = zs21_21 + zvis12 * ze21 + zvis21 * ze12 480 481 ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj) 482 ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj) 483 ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj) 484 ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj) 485 zvis11 = 2.0 * zviseta (ji-1,jj) + dm 486 zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj) 487 zvis12 = zviseta (ji-1,jj) + dm 488 zvis21 = zviseta (ji-1,jj) 489 zdiag = zvis22 * ( ze11 + ze22 ) 490 zs11_12 = zs11_12 + zvis11 * ze11 + zdiag 491 zs12_12 = zs12_12 + zvis12 * ze12 + zvis21 * ze21 492 zs22_12 = zs22_12 + zvis11 * ze22 + zdiag 493 zs21_12 = zs21_12 + zvis12 * ze21 + zvis21 * ze12 494 495 zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 & 496 & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 & 497 & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 & 498 & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 & 499 & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 & 500 & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 & 501 & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 & 502 & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22 503 504 zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 & 505 & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 & 506 & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 & 507 & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 & 508 & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 & 509 & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 & 510 & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 & 511 & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22 512 513 zur = zu_a(ji,jj) - ui_oce(ji,jj) 514 zvr = zv_a(ji,jj) - vi_oce(ji,jj) 515 !!!! 516 zmod = SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 517 za = rhoco * zmod 518 !!!! 519 !!gm chg resul za = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) ) 520 zac = za * cangvg 521 zmpzas = alpha * zcorl(ji,jj) + za * zsang(ji,jj) 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 522 411 zmassdt = zusdtp * zmass(ji,jj) 523 412 zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj) 524 413 525 za1 = zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj) & 526 & + za * ( cangvg * ui_oce(ji,jj) - zsang(ji,jj) * vi_oce(ji,jj) ) 527 za2 = zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj) & 528 & + za * ( cangvg * vi_oce(ji,jj) + zsang(ji,jj) * ui_oce(ji,jj) ) 529 zb1 = zmassdt + zac - zc1u(ji,jj) 530 zb2 = zmpzas - zc2u(ji,jj) 531 zc1 = zmpzas + zc1v(ji,jj) 532 zc2 = zmassdt + zac - zc2v(ji,jj) 533 zdeter = zc1 * zb2 + zc2 * zb1 534 zden = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) ) 535 zunw = ( ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden 536 zvnw = ( ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden 537 zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj) 538 539 zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask 540 zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask 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 ) ) 541 426 END DO 542 427 END DO 543 428 544 CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. ) 545 CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. ) 546 547 ! Test of Convergence 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 501 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 606 END DO 607 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 621 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. 548 627 DO jj = k_j1+1 , k_jpj-1 549 zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) )628 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 550 629 END DO 551 zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) ) 552 !!!! this should be faster on scalar processor 553 ! zresm = MAXVAL( MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ), & 554 ! & ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) ) ) 555 !!!! 630 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 556 631 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 557 632 558 DO jj = k_j1, k_jpj 559 zu_a(:,jj) = zu_n(:,jj) 560 zv_a(:,jj) = zv_n(:,jj) 561 END DO 562 563 IF( zresm <= resl ) EXIT iflag 564 565 ! ! ================ ! 566 END DO iflag ! end Relaxation ! 567 ! ! ================ ! 568 569 IF( zindu == 0 ) THEN ! even iteration 570 DO jj = k_j1 , k_jpj-1 571 zu0(:,jj) = zu_a(:,jj) 572 zv0(:,jj) = zv_a(:,jj) 573 END DO 574 ENDIF 575 ! ! ==================== ! 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 ! ! ==================== ! 576 643 END DO ! end loop over iter ! 577 644 ! ! ==================== ! 578 579 ui_ice(:,:) = zu_a(:,1:jpj)580 vi_ice(:,:) = zv_a(:,1:jpj)581 645 582 646 IF(ln_ctl) THEN 583 647 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter 584 648 CALL prt_ctl_info(charout) 585 CALL prt_ctl(tab2d_1=u i_ice, clinfo1=' lim_rhg : ui_ice :', tab2d_2=vi_ice, clinfo2=' vi_ice :')649 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 586 650 ENDIF 587 651
Note: See TracChangeset
for help on using the changeset viewer.
