limrhg_2.F90 36 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626
  1. MODULE limrhg_2
  2. !!======================================================================
  3. !! *** MODULE limrhg_2 ***
  4. !! Ice rheology : performs sea ice rheology
  5. !!======================================================================
  6. !! History : 0.0 ! 1993-12 (M.A. Morales Maqueda.) Original code
  7. !! 1.0 ! 1994-12 (H. Goosse)
  8. !! 2.0 ! 2003-12 (C. Ethe, G. Madec) F90, mpp
  9. !! - ! 2006-08 (G. Madec) surface module, ice-stress at I-point
  10. !! - ! 2009-09 (G. Madec) Huge verctor optimisation
  11. !! 3.3 ! 2009-05 (G.Garric, C. Bricaud) addition of the lim2_evp case
  12. !!----------------------------------------------------------------------
  13. #if defined key_lim2 && defined key_lim2_vp
  14. !!----------------------------------------------------------------------
  15. !! 'key_lim2' AND LIM-2 sea-ice model
  16. !! 'key_lim2_vp' VP ice rheology
  17. !!----------------------------------------------------------------------
  18. !! lim_rhg_2 : computes ice velocities
  19. !!----------------------------------------------------------------------
  20. USE par_oce ! ocean parameter
  21. USE dom_oce ! ocean space and time domain
  22. USE sbc_oce ! surface boundary condition: ocean variables
  23. USE sbc_ice ! surface boundary condition: ice variables
  24. USE dom_ice_2 ! LIM2: ice domain
  25. USE phycst ! physical constant
  26. USE ice_2 ! LIM2: ice variables
  27. USE lbclnk ! lateral boundary condition - MPP exchanges
  28. USE lib_mpp ! MPP library
  29. USE wrk_nemo ! work arrays
  30. USE in_out_manager ! I/O manager
  31. USE prtctl ! Print control
  32. USE oce , ONLY : snwice_mass, snwice_mass_b
  33. USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
  34. #if defined key_agrif
  35. USE agrif_lim2_interp ! nesting
  36. #endif
  37. IMPLICIT NONE
  38. PRIVATE
  39. PUBLIC lim_rhg_2 ! routine called by lim_dyn
  40. REAL(wp) :: rzero = 0._wp ! constant value: zero
  41. REAL(wp) :: rone = 1._wp ! and one
  42. !! * Substitutions
  43. # include "vectopt_loop_substitute.h90"
  44. !!----------------------------------------------------------------------
  45. !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
  46. !! $Id: limrhg_2.F90 3680 2012-11-27 14:42:24Z rblod $
  47. !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
  48. !!----------------------------------------------------------------------
  49. CONTAINS
  50. SUBROUTINE lim_rhg_2( k_j1, k_jpj )
  51. !!-------------------------------------------------------------------
  52. !! *** SUBROUTINR lim_rhg_2 ***
  53. !!
  54. !! ** purpose : determines the velocity field of sea ice by using
  55. !! atmospheric (wind stress) and oceanic (water stress and surface
  56. !! tilt) forcings. Ice-ice interaction is described by a non-linear
  57. !! viscous-plastic law including shear strength and a bulk rheology.
  58. !!
  59. !! ** Action : - compute u_ice, v_ice the sea-ice velocity defined
  60. !! at I-point
  61. !!-------------------------------------------------------------------
  62. INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation
  63. INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation
  64. !!
  65. INTEGER :: ji, jj ! dummy loop indices
  66. INTEGER :: iter, jter ! temporary integers
  67. CHARACTER (len=50) :: charout
  68. REAL(wp) :: ze11 , ze12 , ze22 , ze21 ! local scalars
  69. REAL(wp) :: zt11 , zt12 , zt21 , zt22 ! - -
  70. REAL(wp) :: zvis11, zvis21, zvis12, zvis22 ! - -
  71. REAL(wp) :: zgphsx, ztagnx, zgsshx, zunw, zur, zusw ! - -
  72. REAL(wp) :: zgphsy, ztagny, zgsshy, zvnw, zvr ! - -
  73. REAL(wp) :: zresm, za, zac, zmod
  74. REAL(wp) :: zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal
  75. REAL(wp) :: ztrace2, zdeter, zdelta, zmask, zdgp, zdgi, zdiag
  76. REAL(wp) :: za1, zb1, zc1, zd1
  77. REAL(wp) :: za2, zb2, zc2, zd2, zden
  78. REAL(wp) :: zs11_11, zs11_12, zs11_21, zs11_22
  79. REAL(wp) :: zs12_11, zs12_12, zs12_21, zs12_22
  80. REAL(wp) :: zs21_11, zs21_12, zs21_21, zs21_22
  81. REAL(wp) :: zs22_11, zs22_12, zs22_21, zs22_22
  82. REAL(wp) :: zintb, zintn
  83. REAL(wp), POINTER, DIMENSION(:,:) :: zfrld, zmass, zcorl
  84. REAL(wp), POINTER, DIMENSION(:,:) :: za1ct, za2ct, zresr
  85. REAL(wp), POINTER, DIMENSION(:,:) :: zc1u, zc1v, zc2u, zc2v
  86. REAL(wp), POINTER, DIMENSION(:,:) :: zsang, zpice
  87. REAL(wp), POINTER, DIMENSION(:,:) :: zu0, zv0
  88. REAL(wp), POINTER, DIMENSION(:,:) :: zu_n, zv_n
  89. REAL(wp), POINTER, DIMENSION(:,:) :: zu_a, zv_a
  90. REAL(wp), POINTER, DIMENSION(:,:) :: zviszeta, zviseta
  91. REAL(wp), POINTER, DIMENSION(:,:) :: zzfrld, zztms
  92. REAL(wp), POINTER, DIMENSION(:,:) :: zi1, zi2, zmasst, zpresh
  93. !!-------------------------------------------------------------------
  94. CALL wrk_alloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr )
  95. CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice )
  96. CALL wrk_alloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 )
  97. CALL wrk_alloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 )
  98. ! Store initial velocities
  99. ! ----------------
  100. zztms(:,0 ) = 0._wp ; zzfrld(:,0 ) = 0._wp
  101. zztms(:,jpj+1) = 0._wp ; zzfrld(:,jpj+1) = 0._wp
  102. zu0 (:,0 ) = 0._wp ; zv0 (:,0 ) = 0._wp
  103. zu0 (:,jpj+1) = 0._wp ; zv0 (:,jpj+1) = 0._wp
  104. zztms(:,1:jpj) = tms (:,:) ; zzfrld(:,1:jpj) = frld (:,:)
  105. zu0 (:,1:jpj) = u_ice(:,:) ; zv0 (:,1:jpj) = v_ice(:,:)
  106. zu_a (:, : ) = zu0 (:,:) ; zv_a (:, : ) = zv0 (:,:)
  107. zu_n (:, : ) = zu0 (:,:) ; zv_n (:, : ) = zv0 (:,:)
  108. !i
  109. zi1 (:,:) = 0._wp
  110. zi2 (:,:) = 0._wp
  111. zpresh(:,:) = 0._wp
  112. zmasst(:,:) = 0._wp
  113. !i
  114. !!gm violant
  115. zfrld(:,:) =0._wp
  116. zcorl(:,:) =0._wp
  117. zmass(:,:) =0._wp
  118. za1ct(:,:) =0._wp
  119. za2ct(:,:) =0._wp
  120. !!gm end
  121. zviszeta(:,:) = 0._wp
  122. zviseta (:,:) = 0._wp
  123. !i zviszeta(:,0 ) = 0._wp ; zviseta(:,0 ) = 0._wp
  124. !i zviszeta(:,jpj ) = 0._wp ; zviseta(:,jpj ) = 0._wp
  125. !i zviszeta(:,jpj+1) = 0._wp ; zviseta(:,jpj+1) = 0._wp
  126. IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
  127. !
  128. ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
  129. ! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
  130. zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
  131. !
  132. ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
  133. ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
  134. zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
  135. !
  136. zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
  137. !
  138. !
  139. ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
  140. zpice(:,:) = ssh_m(:,:)
  141. ENDIF
  142. #if defined key_agrif
  143. ! load the boundary value of velocity in special array zuive and zvice
  144. CALL agrif_rhg_lim2_load
  145. #endif
  146. ! Ice mass, ice strength, and wind stress at the center |
  147. ! of the grid squares. |
  148. !-------------------------------------------------------------------
  149. !CDIR NOVERRCHK
  150. DO jj = k_j1 , k_jpj-1
  151. !CDIR NOVERRCHK
  152. DO ji = 1 , jpi
  153. ! only the sinus changes its sign with the hemisphere
  154. zsang(ji,jj) = SIGN( 1._wp, fcor(ji,jj) ) * sangvg ! only the sinus changes its sign with the hemisphere
  155. !
  156. zmasst(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
  157. zpresh(ji,jj) = tms(ji,jj) * pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
  158. !!gm :: stress given at I-point (F-point for the ocean) only compute the ponderation with the ice fraction (1-frld)
  159. zi1(ji,jj) = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
  160. zi2(ji,jj) = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
  161. END DO
  162. END DO
  163. !---------------------------------------------------------------------------
  164. ! Wind stress, coriolis and mass terms at the corners of the grid squares |
  165. ! Gradient of ice strenght. |
  166. !---------------------------------------------------------------------------
  167. DO jj = k_j1+1, k_jpj-1
  168. DO ji = 2, jpi ! NO vector opt.
  169. zstms = zztms(ji,jj ) * wght(ji,jj,2,2) + zztms(ji-1,jj ) * wght(ji,jj,1,2) &
  170. & + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1)
  171. zusw = 1._wp / MAX( zstms, epsd )
  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)
  176. ! Leads area.
  177. zfrld(ji,jj) = ( zt11 * wght(ji,jj,2,2) + zt12 * wght(ji,jj,1,2) &
  178. & + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw
  179. ! Mass and coriolis coeff. at I-point
  180. zmass(ji,jj) = ( zmasst(ji,jj ) * wght(ji,jj,2,2) + zmasst(ji-1,jj ) * wght(ji,jj,1,2) &
  181. & + zmasst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
  182. zcorl(ji,jj) = zmass(ji,jj) &
  183. & *( fcor(ji,jj ) * wght(ji,jj,2,2) + fcor(ji-1,jj )*wght(ji,jj,1,2) &
  184. & + fcor(ji,jj-1) * wght(ji,jj,2,1) + fcor(ji-1,jj-1)*wght(ji,jj,1,1) ) * zusw
  185. ! Wind stress.
  186. ! always provide stress at I-point
  187. ztagnx = ( zi1(ji,jj ) * wght(ji,jj,2,2) + zi1(ji-1,jj ) * wght(ji,jj,1,2) &
  188. & + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utau_ice(ji,jj)
  189. ztagny = ( zi2(ji,jj ) * wght(ji,jj,2,2) + zi2(ji-1,jj ) * wght(ji,jj,1,2) &
  190. & + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtau_ice(ji,jj)
  191. ! Gradient of ice strength
  192. zgphsx = ( alambd(ji,jj,2,2,2,1) - alambd(ji,jj,2,1,2,1) ) * zpresh(ji ,jj-1) &
  193. & + ( alambd(ji,jj,2,2,2,2) - alambd(ji,jj,2,1,2,2) ) * zpresh(ji ,jj ) &
  194. & - ( alambd(ji,jj,2,2,1,1) + alambd(ji,jj,2,1,1,1) ) * zpresh(ji-1,jj-1) &
  195. & - ( alambd(ji,jj,2,2,1,2) + alambd(ji,jj,2,1,1,2) ) * zpresh(ji-1,jj )
  196. zgphsy = - ( alambd(ji,jj,1,1,2,1) + alambd(ji,jj,1,2,2,1) ) * zpresh(ji ,jj-1) &
  197. & - ( alambd(ji,jj,1,1,1,1) + alambd(ji,jj,1,2,1,1) ) * zpresh(ji-1,jj-1) &
  198. & + ( alambd(ji,jj,1,1,2,2) - alambd(ji,jj,1,2,2,2) ) * zpresh(ji ,jj ) &
  199. & + ( alambd(ji,jj,1,1,1,2) - alambd(ji,jj,1,2,1,2) ) * zpresh(ji-1,jj )
  200. ! Gradient of the sea surface height
  201. zgsshx = ( (zpice(ji ,jj ) - zpice(ji-1,jj ))/e1u(ji-1,jj ) &
  202. & + (zpice(ji ,jj-1) - zpice(ji-1,jj-1))/e1u(ji-1,jj-1) ) * 0.5_wp
  203. zgsshy = ( (zpice(ji ,jj ) - zpice(ji ,jj-1))/e2v(ji ,jj-1) &
  204. & + (zpice(ji-1,jj ) - zpice(ji-1,jj-1))/e2v(ji-1,jj-1) ) * 0.5_wp
  205. ! Computation of the velocity field taking into account the ice-ice interaction.
  206. ! Terms that are independent of the ice velocity field.
  207. za1ct(ji,jj) = ztagnx - zmass(ji,jj) * grav * zgsshx - zgphsx
  208. za2ct(ji,jj) = ztagny - zmass(ji,jj) * grav * zgsshy - zgphsy
  209. END DO
  210. END DO
  211. ! SOLUTION OF THE MOMENTUM EQUATION.
  212. !------------------------------------------
  213. ! ! ==================== !
  214. DO iter = 1 , 2 * nbiter ! loop over iter !
  215. ! ! ==================== !
  216. zindu = MOD( iter , 2 )
  217. zusdtp = ( zindu * 2._wp + ( 1._wp - zindu ) * 1._wp ) * REAL( nbiter ) / rdt_ice
  218. ! Computation of free drift field for free slip boundary conditions.
  219. !CDIR NOVERRCHK
  220. DO jj = k_j1, k_jpj-1
  221. !CDIR NOVERRCHK
  222. DO ji = 1, fs_jpim1
  223. !- Rate of strain tensor.
  224. 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) ) &
  225. & + 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) )
  226. 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) ) &
  227. & - 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) )
  228. 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) ) &
  229. & + 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) )
  230. 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) ) &
  231. & - 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) )
  232. !- Rate of strain tensor.
  233. zdgp = zt11 + zt22
  234. zdgi = zt12 + zt21
  235. ztrace2 = zdgp * zdgp
  236. zdeter = zt11 * zt22 - 0.25_wp * zdgi * zdgi
  237. ! Creep limit depends on the size of the grid.
  238. zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4._wp * zdeter ) * usecc2 ), rn_creepl)
  239. !- Computation of viscosities.
  240. zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn )
  241. zviseta (ji,jj) = zviszeta(ji,jj) * usecc2
  242. END DO
  243. END DO
  244. !- Determination of zc1u, zc2u, zc1v and zc2v.
  245. DO jj = k_j1+1, k_jpj-1
  246. DO ji = 2, fs_jpim1 ! NO vector opt.
  247. !* zc1u , zc2v
  248. zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
  249. zvis12 = zviseta (ji-1,jj-1) + dm
  250. zvis21 = zviseta (ji-1,jj-1)
  251. zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
  252. zdiag = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) )
  253. zs11_11 = zvis11 * akappa(ji-1,jj-1,1,1) + zdiag
  254. zs12_11 = zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2)
  255. zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2)
  256. zs22_11 = zvis11 * akappa(ji-1,jj-1,2,1) + zdiag
  257. zvis11 = 2._wp * 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. zvis11 = 2._wp * zviseta (ji-1,jj) + dm
  267. zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj)
  268. zvis12 = zviseta (ji-1,jj) + dm
  269. zvis21 = zviseta (ji-1,jj)
  270. zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) )
  271. zs11_12 = zvis11 * akappa(ji-1,jj,1,1) + zdiag
  272. zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2)
  273. zs22_12 = zvis11 * akappa(ji-1,jj,2,1) + zdiag
  274. zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2)
  275. zvis11 = 2._wp * zviseta (ji,jj) + dm
  276. zvis22 = zviszeta(ji,jj) - zviseta(ji,jj)
  277. zvis12 = zviseta (ji,jj) + dm
  278. zvis21 = zviseta (ji,jj)
  279. zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) )
  280. zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag
  281. zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2)
  282. zs22_22 = zvis11 * akappa(ji,jj,2,1) + zdiag
  283. zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2)
  284. zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 &
  285. & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 &
  286. & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 &
  287. & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 &
  288. & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 &
  289. & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 &
  290. & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 &
  291. & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
  292. zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 &
  293. & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 &
  294. & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 &
  295. & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 &
  296. & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 &
  297. & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 &
  298. & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 &
  299. & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
  300. !* zc1v , zc2v.
  301. zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
  302. zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
  303. zvis12 = zviseta (ji-1,jj-1) + dm
  304. zvis21 = zviseta (ji-1,jj-1)
  305. zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) )
  306. zs11_11 = zvis11 * akappa(ji-1,jj-1,1,2) + zdiag
  307. zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1)
  308. zs22_11 = zvis11 * akappa(ji-1,jj-1,2,2) + zdiag
  309. zs21_11 = zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1)
  310. zvis11 = 2._wp * zviseta (ji,jj-1) + dm
  311. zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1)
  312. zvis12 = zviseta (ji,jj-1) + dm
  313. zvis21 = zviseta (ji,jj-1)
  314. zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) )
  315. zs11_21 = zvis11 * akappa(ji,jj-1,1,2) + zdiag
  316. zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1)
  317. zs22_21 = zvis11 * akappa(ji,jj-1,2,2) + zdiag
  318. zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1)
  319. zvis11 = 2._wp * zviseta (ji-1,jj) + dm
  320. zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj)
  321. zvis12 = zviseta (ji-1,jj) + dm
  322. zvis21 = zviseta (ji-1,jj)
  323. zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) )
  324. zs11_12 = zvis11 * akappa(ji-1,jj,1,2) + zdiag
  325. zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1)
  326. zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag
  327. zs21_12 = zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1)
  328. zvis11 = 2._wp * zviseta (ji,jj) + dm
  329. zvis22 = zviszeta(ji,jj) - zviseta(ji,jj)
  330. zvis12 = zviseta (ji,jj) + dm
  331. zvis21 = zviseta (ji,jj)
  332. zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) )
  333. zs11_22 = zvis11 * akappa(ji,jj,1,2) + zdiag
  334. zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1)
  335. zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag
  336. zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1)
  337. zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 &
  338. & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 &
  339. & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 &
  340. & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 &
  341. & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 &
  342. & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 &
  343. & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 &
  344. & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
  345. zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 &
  346. & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 &
  347. & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 &
  348. & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 &
  349. & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 &
  350. & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 &
  351. & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 &
  352. & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
  353. END DO
  354. END DO
  355. ! GAUSS-SEIDEL method
  356. ! ! ================ !
  357. iflag: DO jter = 1 , nbitdr ! Relaxation !
  358. ! ! ================ !
  359. !CDIR NOVERRCHK
  360. DO jj = k_j1+1, k_jpj-1
  361. !CDIR NOVERRCHK
  362. DO ji = 2, fs_jpim1 ! NO vector opt.
  363. !
  364. ze11 = akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj)
  365. ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj)
  366. ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj)
  367. ze21 = akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj)
  368. zvis11 = 2._wp * zviseta (ji,jj-1) + dm
  369. zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1)
  370. zvis12 = zviseta (ji,jj-1) + dm
  371. zvis21 = zviseta (ji,jj-1)
  372. zdiag = zvis22 * ( ze11 + ze22 )
  373. zs11_21 = zvis11 * ze11 + zdiag
  374. zs12_21 = zvis12 * ze12 + zvis21 * ze21
  375. zs22_21 = zvis11 * ze22 + zdiag
  376. zs21_21 = zvis12 * ze21 + zvis21 * ze12
  377. ze11 = akappa(ji-1,jj,1,1) * ( zu_a(ji ,jj+1) - zu_a(ji-1,jj+1) ) &
  378. & + akappa(ji-1,jj,1,2) * ( zv_a(ji ,jj+1) + zv_a(ji-1,jj+1) )
  379. ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) ) &
  380. & - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) )
  381. ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji ,jj+1) ) &
  382. & + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji ,jj+1) )
  383. ze21 = akappa(ji-1,jj,1,1) * ( zv_a(ji ,jj+1) - zv_a(ji-1,jj+1) ) &
  384. & - akappa(ji-1,jj,1,2) * ( zu_a(ji ,jj+1) + zu_a(ji-1,jj+1) )
  385. zvis11 = 2._wp * zviseta (ji-1,jj) + dm
  386. zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj)
  387. zvis12 = zviseta (ji-1,jj) + dm
  388. zvis21 = zviseta (ji-1,jj)
  389. zdiag = zvis22 * ( ze11 + ze22 )
  390. zs11_12 = zvis11 * ze11 + zdiag
  391. zs12_12 = zvis12 * ze12 + zvis21 * ze21
  392. zs22_12 = zvis11 * ze22 + zdiag
  393. zs21_12 = zvis12 * ze21 + zvis21 * ze12
  394. ze11 = akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji ,jj+1) ) &
  395. & + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji ,jj+1) )
  396. ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji ,jj+1) - zu_a(ji+1,jj+1) ) &
  397. & - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji ,jj+1) + zv_a(ji+1,jj+1) )
  398. ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji ,jj+1) - zv_a(ji+1,jj+1) ) &
  399. & + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji ,jj+1) + zu_a(ji+1,jj+1) )
  400. ze21 = akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji ,jj+1) ) &
  401. & - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji ,jj+1) )
  402. zvis11 = 2._wp * zviseta (ji,jj) + dm
  403. zvis22 = zviszeta(ji,jj) - zviseta(ji,jj)
  404. zvis12 = zviseta (ji,jj) + dm
  405. zvis21 = zviseta (ji,jj)
  406. zdiag = zvis22 * ( ze11 + ze22 )
  407. zs11_22 = zvis11 * ze11 + zdiag
  408. zs12_22 = zvis12 * ze12 + zvis21 * ze21
  409. zs22_22 = zvis11 * ze22 + zdiag
  410. zs21_22 = zvis12 * ze21 + zvis21 * ze12
  411. ! 2nd part
  412. ze11 = akappa(ji-1,jj-1,1,1) * ( zu_a(ji ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) ) &
  413. & + akappa(ji-1,jj-1,1,2) * ( zv_a(ji ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) )
  414. ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) - zu_a(ji-1,jj) ) &
  415. & - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) + zv_a(ji-1,jj) )
  416. ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji ,jj-1) - zv_a(ji-1,jj) ) &
  417. & + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji ,jj-1) + zu_a(ji-1,jj) )
  418. ze21 = akappa(ji-1,jj-1,1,1) * ( zv_a(ji ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) ) &
  419. & - akappa(ji-1,jj-1,1,2) * ( zu_a(ji ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) )
  420. zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
  421. zvis22 = zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
  422. zvis12 = zviseta (ji-1,jj-1) + dm
  423. zvis21 = zviseta (ji-1,jj-1)
  424. zdiag = zvis22 * ( ze11 + ze22 )
  425. zs11_11 = zvis11 * ze11 + zdiag
  426. zs12_11 = zvis12 * ze12 + zvis21 * ze21
  427. zs22_11 = zvis11 * ze22 + zdiag
  428. zs21_11 = zvis12 * ze21 + zvis21 * ze12
  429. ze11 = akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji ,jj-1) ) &
  430. & + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji ,jj-1) )
  431. ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) ) &
  432. & - akappa(ji,jj-1,2,1) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) )
  433. ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji ,jj-1) + zv_a(ji+1,jj-1) ) &
  434. & + akappa(ji,jj-1,2,1) * ( zu_a(ji ,jj-1) + zu_a(ji+1,jj-1) )
  435. ze21 = akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji ,jj-1) ) &
  436. & - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji ,jj-1) )
  437. zvis11 = 2._wp * zviseta (ji,jj-1) + dm
  438. zvis22 = zviszeta(ji,jj-1) - zviseta(ji,jj-1)
  439. zvis12 = zviseta (ji,jj-1) + dm
  440. zvis21 = zviseta (ji,jj-1)
  441. zdiag = zvis22 * ( ze11 + ze22 )
  442. zs11_21 = zs11_21 + zvis11 * ze11 + zdiag
  443. zs12_21 = zs12_21 + zvis12 * ze12 + zvis21 * ze21
  444. zs22_21 = zs22_21 + zvis11 * ze22 + zdiag
  445. zs21_21 = zs21_21 + zvis12 * ze21 + zvis21 * ze12
  446. ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj)
  447. ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj)
  448. ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj)
  449. ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj)
  450. zvis11 = 2._wp * zviseta (ji-1,jj) + dm
  451. zvis22 = zviszeta(ji-1,jj) - zviseta(ji-1,jj)
  452. zvis12 = zviseta (ji-1,jj) + dm
  453. zvis21 = zviseta (ji-1,jj)
  454. zdiag = zvis22 * ( ze11 + ze22 )
  455. zs11_12 = zs11_12 + zvis11 * ze11 + zdiag
  456. zs12_12 = zs12_12 + zvis12 * ze12 + zvis21 * ze21
  457. zs22_12 = zs22_12 + zvis11 * ze22 + zdiag
  458. zs21_12 = zs21_12 + zvis12 * ze21 + zvis21 * ze12
  459. zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22 &
  460. & - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12 &
  461. & - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11 &
  462. & + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12 &
  463. & + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21 &
  464. & + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22 &
  465. & - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21 &
  466. & - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
  467. zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22 &
  468. & - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12 &
  469. & - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11 &
  470. & + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12 &
  471. & - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21 &
  472. & - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22 &
  473. & + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21 &
  474. & + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
  475. zur = zu_a(ji,jj) - u_oce(ji,jj)
  476. zvr = zv_a(ji,jj) - v_oce(ji,jj)
  477. !!!!
  478. zmod = SQRT( zur*zur + zvr*zvr ) * ( 1._wp - zfrld(ji,jj) )
  479. za = rhoco * zmod
  480. !!!!
  481. !!gm chg resul za = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1._wp - zfrld(ji,jj) )
  482. zac = za * cangvg
  483. zmpzas = alpha * zcorl(ji,jj) + za * zsang(ji,jj)
  484. zmassdt = zusdtp * zmass(ji,jj)
  485. zcorlal = ( 1._wp - alpha ) * zcorl(ji,jj)
  486. za1 = zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj) &
  487. & + za * ( cangvg * u_oce(ji,jj) - zsang(ji,jj) * v_oce(ji,jj) )
  488. za2 = zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj) &
  489. & + za * ( cangvg * v_oce(ji,jj) + zsang(ji,jj) * u_oce(ji,jj) )
  490. zb1 = zmassdt + zac - zc1u(ji,jj)
  491. zb2 = zmpzas - zc2u(ji,jj)
  492. zc1 = zmpzas + zc1v(ji,jj)
  493. zc2 = zmassdt + zac - zc2v(ji,jj)
  494. zdeter = zc1 * zb2 + zc2 * zb1
  495. zden = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) )
  496. zunw = ( ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden
  497. zvnw = ( ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden
  498. zmask = ( 1._wp - MAX( rzero, SIGN( rone , 1._wp - zmass(ji,jj) ) ) ) * tmu(ji,jj)
  499. zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask
  500. zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask
  501. END DO
  502. END DO
  503. CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. )
  504. CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. )
  505. #if defined key_agrif
  506. ! copy the boundary value from u_ice_nst and v_ice_nst to u_ice and v_ice
  507. ! before next interations
  508. CALL agrif_rhg_lim2(zu_n,zv_n)
  509. #endif
  510. ! Test of Convergence
  511. DO jj = k_j1+1 , k_jpj-1
  512. zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) )
  513. END DO
  514. zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) )
  515. !!!! this should be faster on scalar processor
  516. ! zresm = MAXVAL( MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ), &
  517. ! & ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) ) )
  518. !!!!
  519. IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
  520. DO jj = k_j1, k_jpj
  521. zu_a(:,jj) = zu_n(:,jj)
  522. zv_a(:,jj) = zv_n(:,jj)
  523. END DO
  524. IF( zresm <= resl ) EXIT iflag
  525. ! ! ================ !
  526. END DO iflag ! end Relaxation !
  527. ! ! ================ !
  528. IF( zindu == 0 ) THEN ! even iteration
  529. DO jj = k_j1 , k_jpj-1
  530. zu0(:,jj) = zu_a(:,jj)
  531. zv0(:,jj) = zv_a(:,jj)
  532. END DO
  533. ENDIF
  534. ! ! ==================== !
  535. END DO ! end loop over iter !
  536. ! ! ==================== !
  537. u_ice(:,:) = zu_a(:,1:jpj)
  538. v_ice(:,:) = zv_a(:,1:jpj)
  539. IF(ln_ctl) THEN
  540. WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
  541. CALL prt_ctl_info(charout)
  542. CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
  543. ENDIF
  544. CALL wrk_dealloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr )
  545. CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice )
  546. CALL wrk_dealloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 )
  547. CALL wrk_dealloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 )
  548. END SUBROUTINE lim_rhg_2
  549. #else
  550. !!----------------------------------------------------------------------
  551. !! Default option Dummy module NO VP & LIM-2 sea-ice model
  552. !!----------------------------------------------------------------------
  553. CONTAINS
  554. SUBROUTINE lim_rhg_2( k1 , k2 ) ! Dummy routine
  555. WRITE(*,*) 'lim_rhg_2: You should not have seen this print! error?', k1, k2
  556. END SUBROUTINE lim_rhg_2
  557. #endif
  558. !!==============================================================================
  559. END MODULE limrhg_2