New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_opa_interp.F90 on Ticket #107 – Attachment – NEMO

Ticket #107: agrif_opa_interp.F90

File agrif_opa_interp.F90, 61.6 KB (added by nemo_user, 16 years ago)
Line 
1MODULE agrif_opa_interp
2#if defined key_agrif
3   USE par_oce
4   USE oce
5   USE dom_oce     
6   USE sol_oce
7   use dynspg_oce      ! surface pressure gradient variables
8# if defined key_ice_lim
9   USE ice
10   USE ice_oce
11   USE dom_ice
12   USE limtrp
13#endif
14
15   IMPLICIT NONE
16   PRIVATE
17   
18   PUBLIC Agrif_tra, Agrif_dyn, interpu, interpv
19
20# if   defined key_dynspg_ts   ||  defined key_vvl   ||  defined key_esopa
21! FD: added functionality for time-splitting
22   PUBLIC Agrif_dyn_ts
23# endif
24
25# if defined key_ice_lim
26   PUBLIC Agrif_tra_ice, interpu2di, interpv2di
27   PUBLIC Agrif_dyn_ice_load, Agrif_dyn_ice_copy
28#endif
29
30# if defined key_ice_lim
31   REAL(wp), DIMENSION(jpi,jpj), PUBLIC :: zuice, zvice
32#endif
33
34   CONTAINS
35   
36   SUBROUTINE Agrif_tra( kt )
37      !!---------------------------------------------
38      !!   *** ROUTINE Agrif_Tra ***
39      !!---------------------------------------------
40#  include "domzgr_substitute.h90" 
41#  include "vectopt_loop_substitute.h90"
42     
43      INTEGER, INTENT(in) :: kt
44
45      INTEGER :: ji,jj,jk
46      REAL(wp) :: zrhox
47      REAL(wp) :: alpha1, alpha2, alpha3, alpha4
48      REAL(wp) :: alpha5, alpha6, alpha7
49      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zta, zsa
50      !
51      IF(Agrif_Root()) RETURN
52
53      Agrif_SpecialValue=0.
54      Agrif_UseSpecialValue = .TRUE.
55      zta = 0.e0
56      zsa = 0.e0
57
58      CALL Agrif_Bc_variable(zta,tn)
59      CALL Agrif_Bc_variable(zsa,sn)
60      Agrif_UseSpecialValue = .FALSE.
61
62      zrhox = Agrif_Rhox()
63
64      alpha1 = (zrhox-1.)/2.
65      alpha2 = 1.-alpha1
66
67      alpha3 = (zrhox-1)/(zrhox+1)
68      alpha4 = 1.-alpha3
69
70      alpha6 = 2.*(zrhox-1.)/(zrhox+1.)
71      alpha7 = -(zrhox-1)/(zrhox+3)
72      alpha5 = 1. - alpha6 - alpha7
73
74      IF((nbondi == 1).OR.(nbondi == 2)) THEN
75
76         ta(nlci,:,:) = alpha1 * zta(nlci,:,:) + alpha2 * zta(nlci-1,:,:)
77         sa(nlci,:,:) = alpha1 * zsa(nlci,:,:) + alpha2 * zsa(nlci-1,:,:)
78
79         DO jk=1,jpk     
80            DO jj=1,jpj
81               IF (umask(nlci-2,jj,jk).EQ.0.) THEN
82                  ta(nlci-1,jj,jk) = ta(nlci,jj,jk) * tmask(nlci-1,jj,jk)
83                  sa(nlci-1,jj,jk) = sa(nlci,jj,jk) * tmask(nlci-1,jj,jk)
84               ELSE
85                  ta(nlci-1,jj,jk)=(alpha4*ta(nlci,jj,jk)+alpha3*ta(nlci-2,jj,jk))*tmask(nlci-1,jj,jk)
86                  sa(nlci-1,jj,jk)=(alpha4*sa(nlci,jj,jk)+alpha3*sa(nlci-2,jj,jk))*tmask(nlci-1,jj,jk)
87                  IF (un(nlci-2,jj,jk).GT.0.) THEN
88                     ta(nlci-1,jj,jk)=( alpha6*ta(nlci-2,jj,jk)+alpha5*ta(nlci,jj,jk)  &
89                                      + alpha7*ta(nlci-3,jj,jk) ) * tmask(nlci-1,jj,jk)
90                     sa(nlci-1,jj,jk)=( alpha6*sa(nlci-2,jj,jk)+alpha5*sa(nlci,jj,jk)  &
91                                      + alpha7*sa(nlci-3,jj,jk) ) * tmask(nlci-1,jj,jk)
92                  ENDIF
93               ENDIF
94            END DO
95         END DO
96      ENDIF
97
98      IF((nbondj == 1).OR.(nbondj == 2)) THEN
99
100         ta(:,nlcj,:) = alpha1 * zta(:,nlcj,:) + alpha2 * zta(:,nlcj-1,:)
101         sa(:,nlcj,:) = alpha1 * zsa(:,nlcj,:) + alpha2 * zsa(:,nlcj-1,:)
102
103         DO jk=1,jpk     
104            DO ji=1,jpi
105               IF (vmask(ji,nlcj-2,jk).EQ.0.) THEN
106                  ta(ji,nlcj-1,jk) = ta(ji,nlcj,jk) * tmask(ji,nlcj-1,jk)
107                  sa(ji,nlcj-1,jk) = sa(ji,nlcj,jk) * tmask(ji,nlcj-1,jk)
108               ELSE
109                  ta(ji,nlcj-1,jk)=(alpha4*ta(ji,nlcj,jk)+alpha3*ta(ji,nlcj-2,jk))*tmask(ji,nlcj-1,jk)       
110                  sa(ji,nlcj-1,jk)=(alpha4*sa(ji,nlcj,jk)+alpha3*sa(ji,nlcj-2,jk))*tmask(ji,nlcj-1,jk)
111                  IF (vn(ji,nlcj-2,jk) .GT. 0.) THEN
112                     ta(ji,nlcj-1,jk)=( alpha6*ta(ji,nlcj-2,jk)+alpha5*ta(ji,nlcj,jk)  &
113                                      + alpha7*ta(ji,nlcj-3,jk) ) * tmask(ji,nlcj-1,jk)
114                     sa(ji,nlcj-1,jk)=( alpha6*sa(ji,nlcj-2,jk)+alpha5*sa(ji,nlcj,jk)  &
115                                      + alpha7*sa(ji,nlcj-3,jk))*tmask(ji,nlcj-1,jk)
116                  ENDIF
117               ENDIF
118            END DO
119         END DO
120      ENDIF
121
122      IF((nbondi == -1).OR.(nbondi == 2)) THEN
123         ta(1,:,:) = alpha1 * zta(1,:,:) + alpha2 * zta(2,:,:)
124         sa(1,:,:) = alpha1 * zsa(1,:,:) + alpha2 * zsa(2,:,:)     
125         DO jk=1,jpk     
126            DO jj=1,jpj
127               IF (umask(2,jj,jk).EQ.0.) THEN
128                  ta(2,jj,jk) = ta(1,jj,jk) * tmask(2,jj,jk)
129                  sa(2,jj,jk) = sa(1,jj,jk) * tmask(2,jj,jk)
130               ELSE
131                  ta(2,jj,jk)=(alpha4*ta(1,jj,jk)+alpha3*ta(3,jj,jk))*tmask(2,jj,jk)       
132                  sa(2,jj,jk)=(alpha4*sa(1,jj,jk)+alpha3*sa(3,jj,jk))*tmask(2,jj,jk)
133                  IF (un(2,jj,jk).LT.0.) THEN
134                     ta(2,jj,jk)=(alpha6*ta(3,jj,jk)+alpha5*ta(1,jj,jk)+alpha7*ta(4,jj,jk))*tmask(2,jj,jk)
135                     sa(2,jj,jk)=(alpha6*sa(3,jj,jk)+alpha5*sa(1,jj,jk)+alpha7*sa(4,jj,jk))*tmask(2,jj,jk)
136                  ENDIF
137               ENDIF
138            END DO
139         END DO
140      ENDIF
141
142      IF((nbondj == -1).OR.(nbondj == 2)) THEN
143         ta(:,1,:) = alpha1 * zta(:,1,:) + alpha2 * zta(:,2,:)
144         sa(:,1,:) = alpha1 * zsa(:,1,:) + alpha2 * zsa(:,2,:)
145         DO jk=1,jpk     
146            DO ji=1,jpi
147               IF (vmask(ji,2,jk).EQ.0.) THEN
148                  ta(ji,2,jk)=ta(ji,1,jk) * tmask(ji,2,jk)
149                  sa(ji,2,jk)=sa(ji,1,jk) * tmask(ji,2,jk)
150               ELSE
151                  ta(ji,2,jk)=(alpha4*ta(ji,1,jk)+alpha3*ta(ji,3,jk))*tmask(ji,2,jk)
152                  sa(ji,2,jk)=(alpha4*sa(ji,1,jk)+alpha3*sa(ji,3,jk))*tmask(ji,2,jk) 
153                  IF (vn(ji,2,jk) .LT. 0.) THEN
154                     ta(ji,2,jk)=(alpha6*ta(ji,3,jk)+alpha5*ta(ji,1,jk)+alpha7*ta(ji,4,jk))*tmask(ji,2,jk)
155                     sa(ji,2,jk)=(alpha6*sa(ji,3,jk)+alpha5*sa(ji,1,jk)+alpha7*sa(ji,4,jk))*tmask(ji,2,jk)
156                  ENDIF
157               ENDIF
158            END DO
159         END DO
160      ENDIF
161     
162      ta(:,:,:) = ta(:,:,:) * tmask(:,:,:)
163      sa(:,:,:) = sa(:,:,:) * tmask(:,:,:)
164
165   END SUBROUTINE Agrif_tra
166
167   SUBROUTINE Agrif_dyn( kt )
168      !!---------------------------------------------
169      !!   *** ROUTINE Agrif_DYN ***
170      !!---------------------------------------------
171      USE phycst
172      USE in_out_manager
173
174#  include "domzgr_substitute.h90"
175     
176      INTEGER, INTENT(in) :: kt
177
178      REAL(wp) :: timeref
179      REAL(wp) :: z2dt, znugdt
180      REAL(wp) :: zrhox, rhoy
181      REAL(wp), DIMENSION(jpi,jpj) :: zua2d, zva2d
182      REAL(wp), DIMENSION(jpi,jpj) :: spgu1,spgv1
183      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zua, zva
184      INTEGER :: ji,jj,jk
185
186      IF (Agrif_Root()) RETURN
187
188      zrhox = Agrif_Rhox()
189      rhoy = Agrif_Rhoy()
190
191      timeref = 1.
192
193      ! time step: leap-frog
194      z2dt = 2. * rdt
195      ! time step: Euler if restart from rest
196      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
197      ! coefficients
198      znugdt =  rnu * grav * z2dt   
199
200      Agrif_SpecialValue=0.
201      Agrif_UseSpecialValue = .TRUE.
202      zua = 0.
203      zva = 0.
204      CALL Agrif_Bc_variable(zua,un,procname=interpu)
205      CALL Agrif_Bc_variable(zva,vn,procname=interpv)
206      zua2d = 0.
207      zva2d = 0.
208
209      Agrif_SpecialValue=0.
210      Agrif_UseSpecialValue = .TRUE.
211      CALL Agrif_Bc_variable(zua2d,e1u,calledweight=1.,procname=interpu2d)
212      CALL Agrif_Bc_variable(zva2d,e2v,calledweight=1.,procname=interpv2d)
213      Agrif_UseSpecialValue = .FALSE.
214
215
216      IF((nbondi == -1).OR.(nbondi == 2)) THEN
217
218         DO jj=1,jpj
219            laplacu(2,jj) = timeref * (zua2d(2,jj)/(rhoy*e2u(2,jj)))*umask(2,jj,1)
220         END DO
221
222         DO jk=1,jpkm1
223            DO jj=1,jpj
224               ua(1:2,jj,jk) = (zua(1:2,jj,jk)/(rhoy*e2u(1:2,jj)))
225#if ! defined key_zco
226               ua(1:2,jj,jk) = ua(1:2,jj,jk) / fse3u(1:2,jj,jk)
227#endif
228            END DO
229         END DO
230
231         DO jk=1,jpkm1
232            DO jj=1,jpj
233               ua(2,jj,jk) = (ua(2,jj,jk) - z2dt * znugdt * laplacu(2,jj))*umask(2,jj,jk)
234            END DO
235         END DO
236
237         spgu(2,:)=0.
238
239         DO jk=1,jpkm1
240            DO jj=1,jpj
241               spgu(2,jj)=spgu(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk)
242            END DO
243         END DO
244
245         DO jj=1,jpj
246            IF (umask(2,jj,1).NE.0.) THEN
247               spgu(2,jj)=spgu(2,jj)/hu(2,jj)
248            ENDIF
249         END DO
250
251         DO jk=1,jpkm1
252            DO jj=1,jpj
253               ua(2,jj,jk) = 0.25*(ua(1,jj,jk)+2.*ua(2,jj,jk)+ua(3,jj,jk))
254               ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk)
255            END DO
256         END DO
257
258         spgu1(2,:)=0.
259
260         DO jk=1,jpkm1
261            DO jj=1,jpj
262               spgu1(2,jj)=spgu1(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk)
263            END DO
264         END DO
265
266         DO jj=1,jpj
267            IF (umask(2,jj,1).NE.0.) THEN
268               spgu1(2,jj)=spgu1(2,jj)/hu(2,jj)
269            ENDIF
270         END DO
271
272         DO jk=1,jpkm1
273            DO jj=1,jpj
274               ua(2,jj,jk) = (ua(2,jj,jk)+spgu(2,jj)-spgu1(2,jj))*umask(2,jj,jk)
275            END DO
276         END DO
277
278         DO jk=1,jpkm1
279            DO jj=1,jpj
280               va(2,jj,jk) = (zva(2,jj,jk)/(zrhox*e1v(2,jj)))*vmask(2,jj,jk)
281#if ! defined key_zco
282               va(2,jj,jk) = va(2,jj,jk) / fse3v(2,jj,jk)
283#endif           
284            END DO
285         END DO
286
287         sshn(2,:)=sshn(3,:)*tmask(2,:,1)
288         sshb(2,:)=sshb(3,:)*tmask(2,:,1)
289! FD debug add wn:
290         wn(2,:,1)=wn(3,:,1)*tmask(2,:,1)
291
292      ENDIF
293
294      IF((nbondi == 1).OR.(nbondi == 2)) THEN
295
296         DO jj=1,jpj
297            laplacu(nlci-2,jj) = timeref * (zua2d(nlci-2,jj)/(rhoy*e2u(nlci-2,jj)))
298         END DO
299
300         DO jk=1,jpkm1
301            DO jj=1,jpj
302               ua(nlci-2:nlci-1,jj,jk) = (zua(nlci-2:nlci-1,jj,jk)/(rhoy*e2u(nlci-2:nlci-1,jj)))
303
304#if ! defined key_zco
305               ua(nlci-2:nlci-1,jj,jk) = ua(nlci-2:nlci-1,jj,jk) / fse3u(nlci-2:nlci-1,jj,jk)
306#endif
307
308            END DO
309         END DO
310
311         DO jk=1,jpkm1
312            DO jj=1,jpj
313               ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)- z2dt * znugdt * laplacu(nlci-2,jj))*umask(nlci-2,jj,jk)
314            END DO
315         END DO
316
317
318         spgu(nlci-2,:)=0.
319
320         do jk=1,jpkm1
321            do jj=1,jpj
322               spgu(nlci-2,jj)=spgu(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)
323            enddo
324         enddo
325
326         DO jj=1,jpj
327            IF (umask(nlci-2,jj,1).NE.0.) THEN
328               spgu(nlci-2,jj)=spgu(nlci-2,jj)/hu(nlci-2,jj)
329            ENDIF
330         END DO
331
332         DO jk=1,jpkm1
333            DO jj=1,jpj
334               ua(nlci-2,jj,jk) = 0.25*(ua(nlci-3,jj,jk)+2.*ua(nlci-2,jj,jk)+ua(nlci-1,jj,jk))
335
336               ua(nlci-2,jj,jk) = ua(nlci-2,jj,jk) * umask(nlci-2,jj,jk)
337
338            END DO
339         END DO
340
341         spgu1(nlci-2,:)=0.
342
343         DO jk=1,jpkm1
344            DO jj=1,jpj
345               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)*umask(nlci-2,jj,jk)
346            END DO
347         END DO
348
349         DO jj=1,jpj
350            IF (umask(nlci-2,jj,1).NE.0.) THEN
351               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)/hu(nlci-2,jj)
352            ENDIF
353         END DO
354
355         DO jk=1,jpkm1
356            DO jj=1,jpj
357               ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)+spgu(nlci-2,jj)-spgu1(nlci-2,jj))*umask(nlci-2,jj,jk)
358            END DO
359         END DO
360
361         DO jk=1,jpkm1
362            DO jj=1,jpj-1
363               va(nlci-1,jj,jk) = (zva(nlci-1,jj,jk)/(zrhox*e1v(nlci-1,jj)))*vmask(nlci-1,jj,jk)
364#if ! defined key_zco
365               va(nlci-1,jj,jk) = va(nlci-1,jj,jk) / fse3v(nlci-1,jj,jk)
366#endif
367            END DO
368         END DO
369
370         sshn(nlci-1,:)=sshn(nlci-2,:)*tmask(nlci-1,:,1)
371         sshb(nlci-1,:)=sshb(nlci-2,:)*tmask(nlci-1,:,1)
372! FD debug add wn:
373         wn(nlci-1,:,1)=wn(nlci-2,:,1)*tmask(nlci-1,:,1)
374      ENDIF
375
376      IF((nbondj == -1).OR.(nbondj == 2)) THEN
377
378         DO ji=1,jpi
379            laplacv(ji,2) = timeref * (zva2d(ji,2)/(zrhox*e1v(ji,2)))
380         END DO
381
382         DO jk=1,jpkm1
383            DO ji=1,jpi
384               va(ji,1:2,jk) = (zva(ji,1:2,jk)/(zrhox*e1v(ji,1:2)))
385#if ! defined key_zco
386               va(ji,1:2,jk) = va(ji,1:2,jk) / fse3v(ji,1:2,jk)
387#endif
388            END DO
389         END DO
390
391         DO jk=1,jpkm1
392            DO ji=1,jpi
393               va(ji,2,jk) = (va(ji,2,jk) - z2dt * znugdt * laplacv(ji,2))*vmask(ji,2,jk)
394            END DO
395         END DO
396
397         spgv(:,2)=0.
398
399         DO jk=1,jpkm1
400            DO ji=1,jpi
401               spgv(ji,2)=spgv(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)
402            END DO
403         END DO
404
405         DO ji=1,jpi
406            IF (vmask(ji,2,1).NE.0.) THEN
407               spgv(ji,2)=spgv(ji,2)/hv(ji,2)
408            ENDIF
409         END DO
410
411         DO jk=1,jpkm1
412            DO ji=1,jpi
413               va(ji,2,jk)=0.25*(va(ji,1,jk)+2.*va(ji,2,jk)+va(ji,3,jk))
414               va(ji,2,jk)=va(ji,2,jk)*vmask(ji,2,jk)
415            END DO
416         END DO
417
418         spgv1(:,2)=0.
419
420         DO jk=1,jpkm1
421            DO ji=1,jpi
422               spgv1(ji,2)=spgv1(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)*vmask(ji,2,jk)
423            END DO
424         END DO
425
426         DO ji=1,jpi
427            IF (vmask(ji,2,1).NE.0.) THEN
428               spgv1(ji,2)=spgv1(ji,2)/hv(ji,2)
429            ENDIF
430         END DO
431
432         DO jk=1,jpkm1
433            DO ji=1,jpi
434               va(ji,2,jk) = (va(ji,2,jk)+spgv(ji,2)-spgv1(ji,2))*vmask(ji,2,jk)
435            END DO
436         END DO
437
438         DO jk=1,jpkm1
439            DO ji=1,jpi
440               ua(ji,2,jk) = (zua(ji,2,jk)/(rhoy*e2u(ji,2)))*umask(ji,2,jk) 
441#if ! defined key_zco
442               ua(ji,2,jk) = ua(ji,2,jk) / fse3u(ji,2,jk)
443#endif               
444            END DO
445         END DO
446
447         sshn(:,2)=sshn(:,3)*tmask(:,2,1)
448         sshb(:,2)=sshb(:,3)*tmask(:,2,1)
449! FD debug add wn:
450         wn(:,2,1)=wn(:,3,1)*tmask(:,2,1)
451      ENDIF
452
453      IF((nbondj == 1).OR.(nbondj == 2)) THEN
454
455         DO ji=1,jpi
456            laplacv(ji,nlcj-2) = timeref * (zva2d(ji,nlcj-2)/(zrhox*e1v(ji,nlcj-2)))
457         END DO
458
459         DO jk=1,jpkm1
460            DO ji=1,jpi
461               va(ji,nlcj-2:nlcj-1,jk) = (zva(ji,nlcj-2:nlcj-1,jk)/(zrhox*e1v(ji,nlcj-2:nlcj-1)))
462#if ! defined key_zco
463               va(ji,nlcj-2:nlcj-1,jk) = va(ji,nlcj-2:nlcj-1,jk) / fse3v(ji,nlcj-2:nlcj-1,jk)
464#endif
465            END DO
466         END DO
467
468         DO jk=1,jpkm1
469            DO ji=1,jpi
470               va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)-z2dt * znugdt * laplacv(ji,nlcj-2))*vmask(ji,nlcj-2,jk)
471            END DO
472         END DO
473
474
475         spgv(:,nlcj-2)=0.
476
477         DO jk=1,jpkm1
478            DO ji=1,jpi
479               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
480            END DO
481         END DO
482
483         DO ji=1,jpi
484            IF (vmask(ji,nlcj-2,1).NE.0.) THEN
485               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)/hv(ji,nlcj-2)
486            ENDIF
487         END DO
488
489         DO jk=1,jpkm1
490            DO ji=1,jpi
491               va(ji,nlcj-2,jk)=0.25*(va(ji,nlcj-3,jk)+2.*va(ji,nlcj-2,jk)+va(ji,nlcj-1,jk))
492               va(ji,nlcj-2,jk) = va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk)
493            END DO
494         END DO
495
496         spgv1(:,nlcj-2)=0.
497
498         DO jk=1,jpkm1
499            DO ji=1,jpi
500               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
501            END DO
502         END DO
503
504         DO ji=1,jpi
505            IF (vmask(ji,nlcj-2,1).NE.0.) THEN
506               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)/hv(ji,nlcj-2)
507            ENDIF
508         END DO
509
510         DO jk=1,jpkm1
511            DO ji=1,jpi
512               va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)+spgv(ji,nlcj-2)-spgv1(ji,nlcj-2))*vmask(ji,nlcj-2,jk)
513            END DO
514         END DO
515
516         DO jk=1,jpkm1
517            DO ji=1,jpi
518               ua(ji,nlcj-1,jk) = (zua(ji,nlcj-1,jk)/(rhoy*e2u(ji,nlcj-1)))*umask(ji,nlcj-1,jk)
519#if ! defined key_zco
520               ua(ji,nlcj-1,jk) = ua(ji,nlcj-1,jk) / fse3u(ji,nlcj-1,jk)
521#endif         
522            END DO
523         END DO
524
525         sshn(:,nlcj-1)=sshn(:,nlcj-2)*tmask(:,nlcj-1,1)
526         sshb(:,nlcj-1)=sshb(:,nlcj-2)*tmask(:,nlcj-1,1)
527
528! FD debug add wn:
529         wn(:,nlcj-1,1)=wn(:,nlcj-2,1)*tmask(:,nlcj-1,1)
530      ENDIF
531
532   END SUBROUTINE Agrif_dyn
533
534# if   defined key_dynspg_ts   ||  defined key_vvl   ||  defined key_esopa
535!************************************************************************
536! FD: need a special routine for dealing with exchanging data
537!     between the child and parent grid during time-splitting
538!     loop
539!************************************************************************
540   SUBROUTINE Agrif_dyn_ts( kt, kbt, sshb_e )
541      !!---------------------------------------------
542      !!   *** ROUTINE Agrif_DYN ***
543      !!---------------------------------------------
544      USE phycst
545      USE in_out_manager
546
547#  include "domzgr_substitute.h90"
548     
549      INTEGER, INTENT(in) :: kt, kbt
550      REAL(wp), DIMENSION(jpi,jpj) :: sshb_e
551
552      REAL(wp) :: timeref
553      REAL(wp) :: z2dt, znugdt
554      REAL(wp) :: zrhox, rhoy
555      REAL(wp), DIMENSION(jpi,jpj) :: zua2d, zva2d
556      REAL(wp), DIMENSION(jpi,jpj) :: spgu1,spgv1
557      REAL(wp), DIMENSION(jpi,jpj) :: zua, zva
558      INTEGER :: ji,jj
559
560      IF (Agrif_Root()) RETURN
561
562      zrhox = Agrif_Rhox()
563      rhoy  = Agrif_Rhoy()
564
565      timeref = 1.
566
567      ! time step: leap-frog
568      z2dt = 2. * rdtbt
569      IF( kbt == 1 ) z2dt = rdtbt
570      ! coefficients
571      znugdt =  rnu * grav * z2dt   
572
573      Agrif_SpecialValue=0.
574      Agrif_UseSpecialValue = .TRUE.
575      zua = 0.
576      zva = 0.
577      CALL Agrif_Bc_variable(zua,un_b,procname=interpu2d2)
578      CALL Agrif_Bc_variable(zva,vn_b,procname=interpv2d2)
579
580      Agrif_SpecialValue=0.
581      Agrif_UseSpecialValue = .TRUE.
582      zua2d = 0.
583      zva2d = 0.
584      CALL Agrif_Bc_variable(zua2d,e1u,calledweight=1.,procname=interpu2d)
585      CALL Agrif_Bc_variable(zva2d,e2v,calledweight=1.,procname=interpv2d)
586      Agrif_UseSpecialValue = .FALSE.
587
588
589      IF((nbondi == -1).OR.(nbondi == 2)) THEN
590
591         DO jj=1,jpj
592            laplacu(2,jj) = timeref * (zua2d(2,jj)/(rhoy*e2u(2,jj)))*umask(2,jj,1)
593         END DO
594
595            DO jj=1,jpj
596               ua_e(1:2,jj) = (zua(1:2,jj)/(rhoy*e2u(1:2,jj)))
597            END DO
598            DO jj=1,jpj
599               ua_e(2,jj) = (ua_e(2,jj) - z2dt * znugdt * laplacu(2,jj))*umask(2,jj,1)
600            END DO
601
602         spgu(2,:)=0.
603
604            DO jj=1,jpj
605               spgu(2,jj)=spgu(2,jj)+ua_e(2,jj)
606            END DO
607
608            DO jj=1,jpj
609               ua_e(2,jj) = 0.25*(ua_e(1,jj)+2.*ua_e(2,jj)+ua_e(3,jj))
610               ua_e(2,jj) = ua_e(2,jj) * umask(2,jj,1)
611            END DO
612
613         spgu1(2,:)=0.
614
615            DO jj=1,jpj
616               spgu1(2,jj)=spgu1(2,jj)+ua_e(2,jj)
617            END DO
618
619            DO jj=1,jpj
620               ua_e(2,jj) = (ua_e(2,jj)+spgu(2,jj)-spgu1(2,jj))*umask(2,jj,1)
621            END DO
622
623            DO jj=1,jpj
624               va_e(2,jj) = (zva(2,jj)/(zrhox*e1v(2,jj)))*vmask(2,jj,1)
625            END DO
626
627         sshn_e(2,:)=sshn_e(3,:)
628         sshb_e(2,:)=sshb_e(3,:)
629
630      ENDIF
631
632      IF((nbondi == 1).OR.(nbondi == 2)) THEN
633
634         DO jj=1,jpj
635            laplacu(nlci-2,jj) = timeref * (zua2d(nlci-2,jj)/(rhoy*e2u(nlci-2,jj)))
636         END DO
637
638            DO jj=1,jpj
639               ua_e(nlci-2:nlci-1,jj) = (zua(nlci-2:nlci-1,jj)/(rhoy*e2u(nlci-2:nlci-1,jj)))
640            END DO
641            DO jj=1,jpj
642               ua_e(nlci-2,jj) = (ua_e(nlci-2,jj)- z2dt * znugdt * laplacu(nlci-2,jj))*umask(nlci-2,jj,1)
643            END DO
644
645         spgu(nlci-2,:)=0.
646
647            DO jj=1,jpj
648               spgu(nlci-2,jj)=spgu(nlci-2,jj)+ua_e(nlci-2,jj)
649            ENDDO
650
651            DO jj=1,jpj
652               ua_e(nlci-2,jj) = 0.25*(ua_e(nlci-3,jj)+2.*ua_e(nlci-2,jj)+ua_e(nlci-1,jj))
653
654               ua_e(nlci-2,jj) = ua_e(nlci-2,jj) * umask(nlci-2,jj,1)
655
656            END DO
657
658         spgu1(nlci-2,:)=0.
659
660            DO jj=1,jpj
661               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)+ua_e(nlci-2,jj)*umask(nlci-2,jj,1)
662            END DO
663
664            DO jj=1,jpj
665               ua_e(nlci-2,jj) = (ua_e(nlci-2,jj)+spgu(nlci-2,jj)-spgu1(nlci-2,jj))*umask(nlci-2,jj,1)
666            END DO
667            DO jj=1,jpj-1
668               va_e(nlci-1,jj) = (zva(nlci-1,jj)/(zrhox*e1v(nlci-1,jj)))*vmask(nlci-1,jj,1)
669            END DO
670
671         sshn_e(nlci-1,:)=sshn_e(nlci-2,:)
672         sshb_e(nlci-1,:)=sshb_e(nlci-2,:)       
673      ENDIF
674
675      IF((nbondj == -1).OR.(nbondj == 2)) THEN
676
677         DO ji=1,jpi
678            laplacv(ji,2) = timeref * (zva2d(ji,2)/(zrhox*e1v(ji,2)))
679         END DO
680
681            DO ji=1,jpi
682               va_e(ji,1:2) = (zva(ji,1:2)/(zrhox*e1v(ji,1:2)))
683            END DO
684            DO ji=1,jpi
685               va_e(ji,2) = (va_e(ji,2) - z2dt * znugdt * laplacv(ji,2))*vmask(ji,2,1)
686            END DO
687
688         spgv(:,2)=0.
689
690            DO ji=1,jpi
691               spgv(ji,2)=spgv(ji,2)+va_e(ji,2)
692            END DO
693
694            DO ji=1,jpi
695               va_e(ji,2)=0.25*(va_e(ji,1)+2.*va_e(ji,2)+va_e(ji,3))
696               va_e(ji,2)=va_e(ji,2)*vmask(ji,2,1)
697            END DO
698
699         spgv1(:,2)=0.
700
701            DO ji=1,jpi
702               spgv1(ji,2)=spgv1(ji,2)+va_e(ji,2)*vmask(ji,2,1)
703            END DO
704
705            DO ji=1,jpi
706               va_e(ji,2) = (va_e(ji,2)+spgv(ji,2)-spgv1(ji,2))*vmask(ji,2,1)
707            END DO
708            DO ji=1,jpi
709               ua_e(ji,2) = (zua(ji,2)/(rhoy*e2u(ji,2)))*umask(ji,2,1) 
710            END DO
711
712         sshn_e(:,2)=sshn_e(:,3)
713         sshb_e(:,2)=sshb_e(:,3)
714      ENDIF
715
716      IF((nbondj == 1).OR.(nbondj == 2)) THEN
717
718         DO ji=1,jpi
719            laplacv(ji,nlcj-2) = timeref * (zva2d(ji,nlcj-2)/(zrhox*e1v(ji,nlcj-2)))
720         END DO
721
722            DO ji=1,jpi
723               va_e(ji,nlcj-2:nlcj-1) = (zva(ji,nlcj-2:nlcj-1)/(zrhox*e1v(ji,nlcj-2:nlcj-1)))
724            END DO
725            DO ji=1,jpi
726               va_e(ji,nlcj-2) = (va_e(ji,nlcj-2)-z2dt * znugdt * laplacv(ji,nlcj-2))*vmask(ji,nlcj-2,1)
727            END DO
728
729         spgv(:,nlcj-2)=0.
730
731            DO ji=1,jpi
732               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)+va_e(ji,nlcj-2)
733            END DO
734
735            DO ji=1,jpi
736               va_e(ji,nlcj-2)=0.25*(va_e(ji,nlcj-3)+2.*va_e(ji,nlcj-2)+va_e(ji,nlcj-1))
737               va_e(ji,nlcj-2) = va_e(ji,nlcj-2) * vmask(ji,nlcj-2,1)
738            END DO
739
740         spgv1(:,nlcj-2)=0.
741
742            DO ji=1,jpi
743               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)+va_e(ji,nlcj-2)
744            END DO
745
746            DO ji=1,jpi
747               va_e(ji,nlcj-2) = (va_e(ji,nlcj-2)+spgv(ji,nlcj-2)-spgv1(ji,nlcj-2))*vmask(ji,nlcj-2,1)
748            END DO
749            DO ji=1,jpi
750               ua_e(ji,nlcj-1) = (zua(ji,nlcj-1)/(rhoy*e2u(ji,nlcj-1)))*umask(ji,nlcj-1,1)
751            END DO
752
753         sshn_e(:,nlcj-1)=sshn_e(:,nlcj-2)
754         sshb_e(:,nlcj-1)=sshb_e(:,nlcj-2)               
755      ENDIF
756
757   END SUBROUTINE Agrif_dyn_ts
758# endif
759   
760
761# if   defined key_ice_lim
762!************************************************************************
763! FD: need a special routine for dealing with exchanging data
764!     between the child and parent grid during ice step
765!************************************************************************
766   SUBROUTINE Agrif_dyn_ice_load
767      !!---------------------------------------------
768      !!   *** ROUTINE Agrif_DYN ***
769      !!---------------------------------------------
770      USE phycst
771      USE in_out_manager
772
773#  include "domzgr_substitute.h90"
774     
775      IF (Agrif_Root()) RETURN
776
777      Agrif_SpecialValue=0.
778      Agrif_UseSpecialValue = .TRUE.
779      zuice = 0.
780      zvice = 0.
781      CALL Agrif_Bc_variable(zuice,u_ice,procname=interpu2di)
782      CALL Agrif_Bc_variable(zvice,v_ice,procname=interpv2di)
783
784      Agrif_SpecialValue=0.
785      Agrif_UseSpecialValue = .FALSE.
786
787   END SUBROUTINE Agrif_dyn_ice_load
788
789   SUBROUTINE Agrif_dyn_ice_copy
790      !!---------------------------------------------
791      !!   *** ROUTINE Agrif_DYN ***
792      !!---------------------------------------------
793      USE phycst
794      USE in_out_manager
795
796#  include "domzgr_substitute.h90"
797     
798      REAL(wp) :: zrhox, rhoy
799      INTEGER :: ji,jj
800
801      IF (Agrif_Root()) RETURN
802
803      zrhox = Agrif_Rhox()
804      rhoy  = Agrif_Rhoy()
805
806      IF((nbondi == -1).OR.(nbondi == 2)) THEN
807            DO jj=2,jpj
808               u_ice(3,jj) = zuice(3,jj)/(rhoy*e2f(2,jj-1))*tmu(3,jj)
809            END DO
810            DO jj=2,jpj
811               v_ice(3,jj) = zvice(3,jj)/(zrhox*e1f(2,jj-1))*tmu(3,jj)
812            END DO
813      ENDIF
814
815      IF((nbondi == 1).OR.(nbondi == 2)) THEN
816            DO jj=2,jpj
817               u_ice(nlci-1,jj) = zuice(nlci-1,jj)/(rhoy*e2f(nlci-2,jj-1))*tmu(nlci-1,jj)
818            END DO
819            DO jj=2,jpj
820               v_ice(nlci-1,jj) = zvice(nlci-1,jj)/(zrhox*e1f(nlci-2,jj-1))*tmu(nlci-1,jj)
821            END DO
822      ENDIF
823
824      IF((nbondj == -1).OR.(nbondj == 2)) THEN
825            DO ji=2,jpi
826               v_ice(ji,3) = zvice(ji,3)/(zrhox*e1f(ji-1,2))*tmu(ji,3)
827            END DO
828            DO ji=2,jpi
829               u_ice(ji,3) = zuice(ji,3)/(rhoy*e2f(ji-1,2))*tmu(ji,3)
830            END DO
831      ENDIF
832
833      IF((nbondj == 1).OR.(nbondj == 2)) THEN
834            DO ji=2,jpi
835               v_ice(ji,nlcj-1) = zvice(ji,nlcj-1)/(zrhox*e1f(ji-1,nlcj-2))*tmu(ji,nlcj-1)
836            END DO
837            DO ji=2,jpi
838               u_ice(ji,nlcj-1) = zuice(ji,nlcj-1)/(rhoy*e2f(ji-1,nlcj-2))*tmu(ji,nlcj-1)
839            END DO
840      ENDIF
841   END SUBROUTINE Agrif_dyn_ice_copy
842
843   SUBROUTINE Agrif_tra_ice( kt )
844      !!---------------------------------------------
845      !!   *** ROUTINE Agrif_DYN ***
846      !!---------------------------------------------
847      USE phycst
848      USE in_out_manager
849
850#  include "domzgr_substitute.h90"
851     
852      INTEGER, INTENT(in) :: kt
853
854      REAL(wp) :: zrhox, rhoy
855      REAL(wp) :: alpha1, alpha2, alpha3, alpha4
856      REAL(wp) :: alpha5, alpha6, alpha7, zvbord
857      REAL(wp), DIMENSION(jpi,jpj) :: zai, zhi, zhs, zui_u, zvi_v
858      REAL(wp), DIMENSION(jpi,jpj,5) :: zaim, zhim
859      INTEGER :: ji,jj,jk
860
861      IF (Agrif_Root()) RETURN
862
863      zrhox = Agrif_Rhox()
864      rhoy  = Agrif_Rhoy()
865
866      Agrif_SpecialValue=0.
867      Agrif_UseSpecialValue = .TRUE.
868      zai = 0.
869      zhi = 0.
870      zhs = 0.
871
872      CALL Agrif_Bc_variable(zai,frld)
873      CALL Agrif_Bc_variable(zhi,hicif)
874      CALL Agrif_Bc_variable(zhs,hsnif)
875
876! FD add moments for Prather SOM advection
877      zaim=0.
878      zhim=0.
879      CALL Agrif_Bc_variable(zhim(:,:,1),sxice ,procname=interp_sxice)
880      CALL Agrif_Bc_variable(zhim(:,:,2),syice ,procname=interp_syice)
881      CALL Agrif_Bc_variable(zhim(:,:,3),sxxice,procname=interp_sxxice)
882      CALL Agrif_Bc_variable(zhim(:,:,4),syyice,procname=interp_syyice)
883      CALL Agrif_Bc_variable(zhim(:,:,5),sxyice,procname=interp_sxyice)
884      CALL Agrif_Bc_variable(zaim(:,:,1),sxa   ,procname=interp_sxa)
885      CALL Agrif_Bc_variable(zaim(:,:,2),sya   ,procname=interp_sya)
886      CALL Agrif_Bc_variable(zaim(:,:,3),sxxa  ,procname=interp_sxxa)
887      CALL Agrif_Bc_variable(zaim(:,:,4),syya  ,procname=interp_syya)
888      CALL Agrif_Bc_variable(zaim(:,:,5),sxya  ,procname=interp_sxya)
889!FD debug
890if (lwp) write(*,*) 'sxice',sxice(70,jpj-1),sxice(70,jpj-1)/area(70,jpj-1)
891if (lwp) write(*,*) 'zaim ',zaim (70,jpj,1),zaim (70,jpj,1)*area(70,jpj)
892if (lwp) write(*,*) 'zaim ',zaim (70,jpj-1,1),zaim (70,jpj-1,1)*area(70,jpj-1)
893      DO jk=1,5
894        zaim(:,:,jk) = zaim(:,:,jk) * area(:,:)
895      ENDDO
896      DO jk=1,5
897        zhim(:,:,jk) = zhim(:,:,jk) * area(:,:)
898      ENDDO
899
900      Agrif_SpecialValue=0.
901      Agrif_UseSpecialValue = .FALSE.
902
903
904      alpha1 = (zrhox-1.)/2.
905      alpha2 = 1.-alpha1
906
907      alpha3 = (zrhox-1)/(zrhox+1)
908      alpha4 = 1.-alpha3
909
910      alpha6 = 2.*(zrhox-1.)/(zrhox+1.)
911      alpha7 = -(zrhox-1)/(zrhox+3)
912      alpha5 = 1. - alpha6 - alpha7
913
914      zui_u = 0.e0
915      zvi_v = 0.e0
916      ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.       
917      zvbord = 1.0 + ( 1.0 - bound )
918      DO jj = 1, jpjm1
919         DO ji = 1, jpim1
920            zui_u(ji,jj) = ( u_ice(ji+1,jj  ) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) )
921            zvi_v(ji,jj) = ( v_ice(ji  ,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) )
922         END DO
923      END DO
924
925      IF((nbondi == 1).OR.(nbondi == 2)) THEN
926
927         frld (nlci,:) = alpha1 * zai(nlci,:) + alpha2 * zai(nlci-1,:)
928         hicif(nlci,:) = alpha1 * zhi(nlci,:) + alpha2 * zhi(nlci-1,:)
929         hsnif(nlci,:) = alpha1 * zhs(nlci,:) + alpha2 * zhs(nlci-1,:)
930! FD add moments for Prather SOM advection
931         sxice (nlci,:) = alpha1 * zhim(nlci,:,1) + alpha2 * zhim(nlci-1,:,1)
932         syice (nlci,:) = alpha1 * zhim(nlci,:,2) + alpha2 * zhim(nlci-1,:,2)
933         sxxice(nlci,:) = alpha1 * zhim(nlci,:,3) + alpha2 * zhim(nlci-1,:,3)
934         syyice(nlci,:) = alpha1 * zhim(nlci,:,4) + alpha2 * zhim(nlci-1,:,4)
935         sxyice(nlci,:) = alpha1 * zhim(nlci,:,5) + alpha2 * zhim(nlci-1,:,5)
936         sxa   (nlci,:) = alpha1 * zaim(nlci,:,1) + alpha2 * zaim(nlci-1,:,1)
937         sya   (nlci,:) = alpha1 * zaim(nlci,:,2) + alpha2 * zaim(nlci-1,:,2)
938         sxxa  (nlci,:) = alpha1 * zaim(nlci,:,3) + alpha2 * zaim(nlci-1,:,3)
939         syya  (nlci,:) = alpha1 * zaim(nlci,:,4) + alpha2 * zaim(nlci-1,:,4)
940         sxya  (nlci,:) = alpha1 * zaim(nlci,:,5) + alpha2 * zaim(nlci-1,:,5)
941
942         DO jj=1,jpj
943            IF (umask(nlci-2,jj,1).EQ.0.) THEN
944               frld (nlci-1,jj) = frld (nlci,jj) * tms(nlci-1,jj)
945               hicif(nlci-1,jj) = hicif(nlci,jj) * tms(nlci-1,jj)
946               hsnif(nlci-1,jj) = hsnif(nlci,jj) * tms(nlci-1,jj)
947! FD add moments for Prather SOM advection
948               sxice (nlci-1,jj) = sxice (nlci,jj) * tms(nlci-1,jj)
949               syice (nlci-1,jj) = syice (nlci,jj) * tms(nlci-1,jj)
950               sxxice(nlci-1,jj) = sxxice(nlci,jj) * tms(nlci-1,jj)
951               syyice(nlci-1,jj) = syyice(nlci,jj) * tms(nlci-1,jj)
952               sxyice(nlci-1,jj) = sxyice(nlci,jj) * tms(nlci-1,jj)
953               sxa   (nlci-1,jj) = sxa   (nlci,jj) * tms(nlci-1,jj)
954               sya   (nlci-1,jj) = sya   (nlci,jj) * tms(nlci-1,jj)
955               sxxa  (nlci-1,jj) = sxxa  (nlci,jj) * tms(nlci-1,jj)
956               syya  (nlci-1,jj) = syya  (nlci,jj) * tms(nlci-1,jj)
957               sxya  (nlci-1,jj) = sxya  (nlci,jj) * tms(nlci-1,jj)
958
959            ELSE
960               frld (nlci-1,jj)=(alpha4*frld (nlci,jj)+alpha3*frld (nlci-2,jj))*tms(nlci-1,jj)
961               hicif(nlci-1,jj)=(alpha4*hicif(nlci,jj)+alpha3*hicif(nlci-2,jj))*tms(nlci-1,jj)
962               hsnif(nlci-1,jj)=(alpha4*hsnif(nlci,jj)+alpha3*hsnif(nlci-2,jj))*tms(nlci-1,jj)
963! FD add moments for Prather SOM advection
964               sxice (nlci-1,jj)=(alpha4*sxice (nlci,jj)+alpha3*sxice (nlci-2,jj))*tms(nlci-1,jj)
965               syice (nlci-1,jj)=(alpha4*syice (nlci,jj)+alpha3*syice (nlci-2,jj))*tms(nlci-1,jj)
966               sxxice(nlci-1,jj)=(alpha4*sxxice(nlci,jj)+alpha3*sxxice(nlci-2,jj))*tms(nlci-1,jj)
967               syyice(nlci-1,jj)=(alpha4*syyice(nlci,jj)+alpha3*syyice(nlci-2,jj))*tms(nlci-1,jj)
968               sxyice(nlci-1,jj)=(alpha4*sxyice(nlci,jj)+alpha3*sxyice(nlci-2,jj))*tms(nlci-1,jj)
969               sxa   (nlci-1,jj)=(alpha4*sxa   (nlci,jj)+alpha3*sxa   (nlci-2,jj))*tms(nlci-1,jj)
970               sya   (nlci-1,jj)=(alpha4*sya   (nlci,jj)+alpha3*sya   (nlci-2,jj))*tms(nlci-1,jj)
971               sxxa  (nlci-1,jj)=(alpha4*sxxa  (nlci,jj)+alpha3*sxxa  (nlci-2,jj))*tms(nlci-1,jj)
972               syya  (nlci-1,jj)=(alpha4*syya  (nlci,jj)+alpha3*syya  (nlci-2,jj))*tms(nlci-1,jj)
973               sxya  (nlci-1,jj)=(alpha4*sxya  (nlci,jj)+alpha3*sxya  (nlci-2,jj))*tms(nlci-1,jj)
974
975               IF (zui_u(nlci-2,jj).GT.0.) THEN
976              frld (nlci-1,jj)=( alpha6*frld (nlci-2,jj)+alpha5*frld (nlci,jj)  &
977                     + alpha7*frld (nlci-3,jj) ) * tms(nlci-1,jj)
978              hicif(nlci-1,jj)=( alpha6*hicif(nlci-2,jj)+alpha5*hicif(nlci,jj)  &
979                     + alpha7*hicif(nlci-3,jj) ) * tms(nlci-1,jj)
980              hsnif(nlci-1,jj)=( alpha6*hsnif(nlci-2,jj)+alpha5*hsnif(nlci,jj)  &
981                     + alpha7*hsnif(nlci-3,jj) ) * tms(nlci-1,jj)
982! FD add moments for Prather SOM advection
983              sxice (nlci-1,jj)=( alpha6*sxice (nlci-2,jj)+alpha5*sxice (nlci,jj)  &
984                      + alpha7*sxice (nlci-3,jj) ) * tms(nlci-1,jj)
985              syice (nlci-1,jj)=( alpha6*syice (nlci-2,jj)+alpha5*syice (nlci,jj)  &
986                      + alpha7*syice (nlci-3,jj) ) * tms(nlci-1,jj)
987              sxxice(nlci-1,jj)=( alpha6*sxxice(nlci-2,jj)+alpha5*sxxice(nlci,jj)  &
988                           + alpha7*sxxice(nlci-3,jj) ) * tms(nlci-1,jj)
989              syyice(nlci-1,jj)=( alpha6*syyice(nlci-2,jj)+alpha5*syyice(nlci,jj)  &
990                           + alpha7*syyice(nlci-3,jj) ) * tms(nlci-1,jj)
991              sxyice(nlci-1,jj)=( alpha6*sxyice(nlci-2,jj)+alpha5*sxyice(nlci,jj)  &
992                      + alpha7*sxyice(nlci-3,jj) ) * tms(nlci-1,jj)
993              sxa   (nlci-1,jj)=( alpha6*sxa   (nlci-2,jj)+alpha5*sxa   (nlci,jj)  &
994                      + alpha7*sxa   (nlci-3,jj) ) * tms(nlci-1,jj)
995              sya   (nlci-1,jj)=( alpha6*sya   (nlci-2,jj)+alpha5*sya   (nlci,jj)  &
996                      + alpha7*sya   (nlci-3,jj) ) * tms(nlci-1,jj)
997              sxxa  (nlci-1,jj)=( alpha6*sxxa  (nlci-2,jj)+alpha5*sxxa  (nlci,jj)  &
998                           + alpha7*sxxa  (nlci-3,jj) ) * tms(nlci-1,jj)
999              syya  (nlci-1,jj)=( alpha6*syya  (nlci-2,jj)+alpha5*syya  (nlci,jj)  &
1000                           + alpha7*syya  (nlci-3,jj) ) * tms(nlci-1,jj)
1001              sxya  (nlci-1,jj)=( alpha6*sxya  (nlci-2,jj)+alpha5*sxya  (nlci,jj)  &
1002                      + alpha7*sxya  (nlci-3,jj) ) * tms(nlci-1,jj)
1003               ENDIF
1004            ENDIF
1005         END DO
1006!         sxice (nlci-1,:) = sxice (nlci-2,:) * tms(nlci-1,:)
1007!         sxxice(nlci-1,:) = sxxice(nlci-2,:) * tms(nlci-1,:)
1008!         syice (nlci-1,:) = syice (nlci-2,:) * tms(nlci-1,:)
1009!         syyice(nlci-1,:) = syyice(nlci-2,:) * tms(nlci-1,:)
1010!         sxyice(nlci-1,:) = sxyice(nlci-2,:) * tms(nlci-1,:)
1011!         sxsn  (nlci-1,:) = sxsn  (nlci-2,:) * tms(nlci-1,:)
1012!         sxxsn (nlci-1,:) = sxxsn (nlci-2,:) * tms(nlci-1,:)
1013!         sysn  (nlci-1,:) = sysn  (nlci-2,:) * tms(nlci-1,:)
1014!         syysn (nlci-1,:) = syysn (nlci-2,:) * tms(nlci-1,:)
1015!         sxysn (nlci-1,:) = sxysn (nlci-2,:) * tms(nlci-1,:)
1016!         sxa   (nlci-1,:) = sxa   (nlci-2,:) * tms(nlci-1,:)
1017!         sxxa  (nlci-1,:) = sxxa  (nlci-2,:) * tms(nlci-1,:)
1018!         sya   (nlci-1,:) = sya   (nlci-2,:) * tms(nlci-1,:)
1019!         syya  (nlci-1,:) = syya  (nlci-2,:) * tms(nlci-1,:)
1020!         sxya  (nlci-1,:) = sxya  (nlci-2,:) * tms(nlci-1,:)
1021      ENDIF
1022
1023      IF((nbondj == 1).OR.(nbondj == 2)) THEN
1024
1025         frld (:,nlcj) = alpha1 * zai(:,nlcj) + alpha2 * zai(:,nlcj-1)
1026         hicif(:,nlcj) = alpha1 * zhi(:,nlcj) + alpha2 * zhi(:,nlcj-1)
1027         hsnif(:,nlcj) = alpha1 * zhs(:,nlcj) + alpha2 * zhs(:,nlcj-1)
1028! FD add moments for Prather SOM advection
1029         sxice (:,nlcj) = alpha1 * zhim(:,nlcj,1) + alpha2 * zhim(:,nlcj-1,1)
1030         syice (:,nlcj) = alpha1 * zhim(:,nlcj,2) + alpha2 * zhim(:,nlcj-1,2)
1031         sxxice(:,nlcj) = alpha1 * zhim(:,nlcj,3) + alpha2 * zhim(:,nlcj-1,3)
1032         syyice(:,nlcj) = alpha1 * zhim(:,nlcj,4) + alpha2 * zhim(:,nlcj-1,4)
1033         sxyice(:,nlcj) = alpha1 * zhim(:,nlcj,5) + alpha2 * zhim(:,nlcj-1,5)
1034         sxa   (:,nlcj) = alpha1 * zaim(:,nlcj,1) + alpha2 * zaim(:,nlcj-1,1)
1035         sya   (:,nlcj) = alpha1 * zaim(:,nlcj,2) + alpha2 * zaim(:,nlcj-1,2)
1036         sxxa  (:,nlcj) = alpha1 * zaim(:,nlcj,3) + alpha2 * zaim(:,nlcj-1,3)
1037         syya  (:,nlcj) = alpha1 * zaim(:,nlcj,4) + alpha2 * zaim(:,nlcj-1,4)
1038         sxya  (:,nlcj) = alpha1 * zaim(:,nlcj,5) + alpha2 * zaim(:,nlcj-1,5)
1039
1040         DO ji=1,jpi
1041            IF (vmask(ji,nlcj-2,1).EQ.0.) THEN
1042               frld (ji,nlcj-1) = frld (ji,nlcj) * tms(ji,nlcj-1)
1043               hicif(ji,nlcj-1) = hicif(ji,nlcj) * tms(ji,nlcj-1)
1044               hsnif(ji,nlcj-1) = hsnif(ji,nlcj) * tms(ji,nlcj-1)
1045! FD add moments for Prather SOM advection
1046               sxice (ji,nlcj-1) = sxice (ji,nlcj) * tms(ji,nlcj-1)
1047               syice (ji,nlcj-1) = syice (ji,nlcj) * tms(ji,nlcj-1)
1048               sxxice(ji,nlcj-1) = sxxice(ji,nlcj) * tms(ji,nlcj-1)
1049               syyice(ji,nlcj-1) = syyice(ji,nlcj) * tms(ji,nlcj-1)
1050               sxyice(ji,nlcj-1) = sxyice(ji,nlcj) * tms(ji,nlcj-1)
1051               sxa   (ji,nlcj-1) = sxa   (ji,nlcj) * tms(ji,nlcj-1)
1052               sya   (ji,nlcj-1) = sya   (ji,nlcj) * tms(ji,nlcj-1)
1053               sxxa  (ji,nlcj-1) = sxxa  (ji,nlcj) * tms(ji,nlcj-1)
1054               syya  (ji,nlcj-1) = syya  (ji,nlcj) * tms(ji,nlcj-1)
1055               sxya  (ji,nlcj-1) = sxya  (ji,nlcj) * tms(ji,nlcj-1)
1056
1057            ELSE
1058               frld (ji,nlcj-1)=(alpha4*frld (ji,nlcj)+alpha3*frld (ji,nlcj-2))*tms(ji,nlcj-1)     
1059               hicif(ji,nlcj-1)=(alpha4*hicif(ji,nlcj)+alpha3*hicif(ji,nlcj-2))*tms(ji,nlcj-1)     
1060               hsnif(ji,nlcj-1)=(alpha4*hsnif(ji,nlcj)+alpha3*hsnif(ji,nlcj-2))*tms(ji,nlcj-1)     
1061! FD add moments for Prather SOM advection
1062               sxice (ji,nlcj-1)=(alpha4*sxice (ji,nlcj)+alpha3*sxice (ji,nlcj-2))*tms(ji,nlcj-1)     
1063               syice (ji,nlcj-1)=(alpha4*syice (ji,nlcj)+alpha3*syice (ji,nlcj-2))*tms(ji,nlcj-1)     
1064               sxxice(ji,nlcj-1)=(alpha4*sxxice(ji,nlcj)+alpha3*sxxice(ji,nlcj-2))*tms(ji,nlcj-1) 
1065               syyice(ji,nlcj-1)=(alpha4*syyice(ji,nlcj)+alpha3*syyice(ji,nlcj-2))*tms(ji,nlcj-1) 
1066               sxyice(ji,nlcj-1)=(alpha4*sxyice(ji,nlcj)+alpha3*sxyice(ji,nlcj-2))*tms(ji,nlcj-1) 
1067               sxa   (ji,nlcj-1)=(alpha4*sxa   (ji,nlcj)+alpha3*sxa   (ji,nlcj-2))*tms(ji,nlcj-1)     
1068               sya   (ji,nlcj-1)=(alpha4*sya   (ji,nlcj)+alpha3*sya   (ji,nlcj-2))*tms(ji,nlcj-1)     
1069               sxxa  (ji,nlcj-1)=(alpha4*sxxa  (ji,nlcj)+alpha3*sxxa  (ji,nlcj-2))*tms(ji,nlcj-1) 
1070               syya  (ji,nlcj-1)=(alpha4*syya  (ji,nlcj)+alpha3*syya  (ji,nlcj-2))*tms(ji,nlcj-1) 
1071               sxya  (ji,nlcj-1)=(alpha4*sxya  (ji,nlcj)+alpha3*sxya  (ji,nlcj-2))*tms(ji,nlcj-1) 
1072
1073               IF (zvi_v(ji,nlcj-2) .GT. 0.) THEN
1074              frld (ji,nlcj-1)=( alpha6*frld (ji,nlcj-2)+alpha5*frld (ji,nlcj)  &
1075                     + alpha7*frld (ji,nlcj-3) ) * tms(ji,nlcj-1)
1076              hicif(ji,nlcj-1)=( alpha6*hicif(ji,nlcj-2)+alpha5*hicif(ji,nlcj)  &
1077                     + alpha7*hicif(ji,nlcj-3) ) * tms(ji,nlcj-1)
1078              hsnif(ji,nlcj-1)=( alpha6*hsnif(ji,nlcj-2)+alpha5*hsnif(ji,nlcj)  &
1079                     + alpha7*hsnif(ji,nlcj-3) ) * tms(ji,nlcj-1)
1080! FD add moments for Prather SOM advection
1081              sxice (ji,nlcj-1)=( alpha6*sxice (ji,nlcj-2)+alpha5*sxice (ji,nlcj)  &
1082                      + alpha7*sxice (ji,nlcj-3) ) * tms(ji,nlcj-1)
1083              syice (ji,nlcj-1)=( alpha6*syice (ji,nlcj-2)+alpha5*syice (ji,nlcj)  &
1084                      + alpha7*syice (ji,nlcj-3) ) * tms(ji,nlcj-1)
1085              sxxice(ji,nlcj-1)=( alpha6*sxxice(ji,nlcj-2)+alpha5*sxxice(ji,nlcj)  &
1086                      + alpha7*sxxice(ji,nlcj-3) ) * tms(ji,nlcj-1)
1087              syyice(ji,nlcj-1)=( alpha6*syyice(ji,nlcj-2)+alpha5*syyice(ji,nlcj)  &
1088                      + alpha7*syyice(ji,nlcj-3) ) * tms(ji,nlcj-1)
1089              sxyice(ji,nlcj-1)=( alpha6*sxyice(ji,nlcj-2)+alpha5*sxyice(ji,nlcj)  &
1090                      + alpha7*sxyice(ji,nlcj-3) ) * tms(ji,nlcj-1)
1091              sxa   (ji,nlcj-1)=( alpha6*sxa   (ji,nlcj-2)+alpha5*sxa   (ji,nlcj)  &
1092                      + alpha7*sxa   (ji,nlcj-3) ) * tms(ji,nlcj-1)
1093              sya   (ji,nlcj-1)=( alpha6*sya   (ji,nlcj-2)+alpha5*sya   (ji,nlcj)  &
1094                      + alpha7*sya   (ji,nlcj-3) ) * tms(ji,nlcj-1)
1095              sxxa  (ji,nlcj-1)=( alpha6*sxxa  (ji,nlcj-2)+alpha5*sxxa  (ji,nlcj)  &
1096                      + alpha7*sxxa  (ji,nlcj-3) ) * tms(ji,nlcj-1)
1097              syya  (ji,nlcj-1)=( alpha6*syya  (ji,nlcj-2)+alpha5*syya  (ji,nlcj)  &
1098                      + alpha7*syya  (ji,nlcj-3) ) * tms(ji,nlcj-1)
1099              sxya  (ji,nlcj-1)=( alpha6*sxya  (ji,nlcj-2)+alpha5*sxya  (ji,nlcj)  &
1100                      + alpha7*sxya  (ji,nlcj-3) ) * tms(ji,nlcj-1)
1101               ENDIF
1102            ENDIF
1103         END DO
1104         DO ji=1,jpi
1105!           sxice (ji,nlcj-2:nlcj-1) = sxice (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1106!           sxxice(ji,nlcj-2:nlcj-1) = sxxice(ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1107!           syice (ji,nlcj-2:nlcj-1) = syice (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1108!           syyice(ji,nlcj-2:nlcj-1) = syyice(ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1109!           sxyice(ji,nlcj-2:nlcj-1) = sxyice(ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1110!           sxsn  (ji,nlcj-2:nlcj-1) = sxsn  (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1111!           sxxsn (ji,nlcj-2:nlcj-1) = sxxsn (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1112!           sysn  (ji,nlcj-2:nlcj-1) = sysn  (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1113!           syysn (ji,nlcj-2:nlcj-1) = syysn (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1114!           sxysn (ji,nlcj-2:nlcj-1) = sxysn (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1115!           sxa   (ji,nlcj-2:nlcj-1) = sxa   (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1116!           sxxa  (ji,nlcj-2:nlcj-1) = sxxa  (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1117!           sya   (ji,nlcj-2:nlcj-1) = sya   (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1118!           syya  (ji,nlcj-2:nlcj-1) = syya  (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1119!           sxya  (ji,nlcj-2:nlcj-1) = sxya  (ji,nlcj-3) * tms(ji,nlcj-2:nlcj-1)
1120         END DO
1121      ENDIF
1122
1123      IF((nbondi == -1).OR.(nbondi == 2)) THEN
1124         frld (1,:) = alpha1 * zai(1,:) + alpha2 * zai(2,:)
1125         hicif(1,:) = alpha1 * zhi(1,:) + alpha2 * zhi(2,:)
1126         hsnif(1,:) = alpha1 * zhs(1,:) + alpha2 * zhs(2,:)
1127! FD add moments for Prather SOM advection
1128         sxice (1,:) = alpha1 * zhim(1,:,1) + alpha2 * zhim(2,:,1)
1129         syice (1,:) = alpha1 * zhim(1,:,2) + alpha2 * zhim(2,:,2)
1130         sxxice(1,:) = alpha1 * zhim(1,:,3) + alpha2 * zhim(2,:,3)
1131         syyice(1,:) = alpha1 * zhim(1,:,4) + alpha2 * zhim(2,:,4)
1132         sxyice(1,:) = alpha1 * zhim(1,:,5) + alpha2 * zhim(2,:,5)
1133         sxa   (1,:) = alpha1 * zaim(1,:,1) + alpha2 * zaim(2,:,1)
1134         sya   (1,:) = alpha1 * zaim(1,:,2) + alpha2 * zaim(2,:,2)
1135         sxxa  (1,:) = alpha1 * zaim(1,:,3) + alpha2 * zaim(2,:,3)
1136         syya  (1,:) = alpha1 * zaim(1,:,4) + alpha2 * zaim(2,:,4)
1137         sxya  (1,:) = alpha1 * zaim(1,:,5) + alpha2 * zaim(2,:,5)
1138
1139         DO jj=1,jpj
1140            IF (umask(2,jj,1).EQ.0.) THEN
1141               frld (2,jj) = frld (1,jj) * tms(2,jj)
1142               hicif(2,jj) = hicif(1,jj) * tms(2,jj)
1143               hsnif(2,jj) = hsnif(1,jj) * tms(2,jj)
1144! FD add moments for Prather SOM advection
1145               sxice (2,jj) = sxice (1,jj) * tms(2,jj)
1146               syice (2,jj) = syice (1,jj) * tms(2,jj)
1147               sxxice(2,jj) = sxxice(1,jj) * tms(2,jj)
1148               syyice(2,jj) = syyice(1,jj) * tms(2,jj)
1149               sxyice(2,jj) = sxyice(1,jj) * tms(2,jj)
1150               sxa   (2,jj) = sxa   (1,jj) * tms(2,jj)
1151               sya   (2,jj) = sya   (1,jj) * tms(2,jj)
1152               sxxa  (2,jj) = sxxa  (1,jj) * tms(2,jj)
1153               syya  (2,jj) = syya  (1,jj) * tms(2,jj)
1154               sxya  (2,jj) = sxya  (1,jj) * tms(2,jj)
1155
1156            ELSE
1157               frld (2,jj)=(alpha4*frld (1,jj)+alpha3*frld (3,jj))*tms(2,jj)   
1158               hicif(2,jj)=(alpha4*hicif(1,jj)+alpha3*hicif(3,jj))*tms(2,jj)   
1159               hsnif(2,jj)=(alpha4*hsnif(1,jj)+alpha3*hsnif(3,jj))*tms(2,jj)   
1160! FD add moments for Prather SOM advection
1161               sxice (2,jj)=(alpha4*sxice (1,jj)+alpha3*sxice (3,jj))*tms(2,jj)     
1162               syice (2,jj)=(alpha4*syice (1,jj)+alpha3*syice (3,jj))*tms(2,jj)     
1163               sxxice(2,jj)=(alpha4*sxxice(1,jj)+alpha3*sxxice(3,jj))*tms(2,jj) 
1164               syyice(2,jj)=(alpha4*syyice(1,jj)+alpha3*syyice(3,jj))*tms(2,jj) 
1165               sxyice(2,jj)=(alpha4*sxyice(1,jj)+alpha3*sxyice(3,jj))*tms(2,jj) 
1166               sxa   (2,jj)=(alpha4*sxa   (1,jj)+alpha3*sxa   (3,jj))*tms(2,jj)     
1167               sya   (2,jj)=(alpha4*sya   (1,jj)+alpha3*sya   (3,jj))*tms(2,jj)     
1168               sxxa  (2,jj)=(alpha4*sxxa  (1,jj)+alpha3*sxxa  (3,jj))*tms(2,jj) 
1169               syya  (2,jj)=(alpha4*syya  (1,jj)+alpha3*syya  (3,jj))*tms(2,jj) 
1170               sxya  (2,jj)=(alpha4*sxya  (1,jj)+alpha3*sxya  (3,jj))*tms(2,jj) 
1171
1172               IF (zui_u(2,jj).LT.0.) THEN
1173              frld (2,jj)=(alpha6*frld (3,jj)+alpha5*frld (1,jj)+alpha7*frld (4,jj))*tms(2,jj)
1174              hicif(2,jj)=(alpha6*hicif(3,jj)+alpha5*hicif(1,jj)+alpha7*hicif(4,jj))*tms(2,jj)
1175              hsnif(2,jj)=(alpha6*hsnif(3,jj)+alpha5*hsnif(1,jj)+alpha7*hsnif(4,jj))*tms(2,jj)
1176! FD add moments for Prather SOM advection
1177              sxice (2,jj)=(alpha6*sxice (3,jj)+alpha5*sxice (1,jj)+alpha7*sxice (4,jj))*tms(2,jj)
1178              syice (2,jj)=(alpha6*syice (3,jj)+alpha5*syice (1,jj)+alpha7*syice (4,jj))*tms(2,jj)
1179              sxxice(2,jj)=(alpha6*sxxice(3,jj)+alpha5*sxxice(1,jj)+alpha7*sxxice(4,jj))*tms(2,jj)
1180              syyice(2,jj)=(alpha6*syyice(3,jj)+alpha5*syyice(1,jj)+alpha7*syyice(4,jj))*tms(2,jj)
1181              sxyice(2,jj)=(alpha6*sxyice(3,jj)+alpha5*sxyice(1,jj)+alpha7*sxyice(4,jj))*tms(2,jj)
1182              sxa   (2,jj)=(alpha6*sxa   (3,jj)+alpha5*sxa   (1,jj)+alpha7*sxa   (4,jj))*tms(2,jj)
1183              sya   (2,jj)=(alpha6*sya   (3,jj)+alpha5*sya   (1,jj)+alpha7*sya   (4,jj))*tms(2,jj)
1184              sxxa  (2,jj)=(alpha6*sxxa  (3,jj)+alpha5*sxxa  (1,jj)+alpha7*sxxa  (4,jj))*tms(2,jj)
1185              syya  (2,jj)=(alpha6*syya  (3,jj)+alpha5*syya  (1,jj)+alpha7*syya  (4,jj))*tms(2,jj)
1186              sxya  (2,jj)=(alpha6*sxya  (3,jj)+alpha5*sxya  (1,jj)+alpha7*sxya  (4,jj))*tms(2,jj)
1187               ENDIF
1188            ENDIF
1189         END DO
1190!         sxice (2,:) = sxice (3,:) * tms(2,:)
1191!         sxxice(2,:) = sxxice(3,:) * tms(2,:)
1192!         syice (2,:) = syice (3,:) * tms(2,:)
1193!         syyice(2,:) = syyice(3,:) * tms(2,:)
1194!         sxyice(2,:) = sxyice(3,:) * tms(2,:)
1195!         sxsn  (2,:) = sxsn  (3,:) * tms(2,:)
1196!         sxxsn (2,:) = sxxsn (3,:) * tms(2,:)
1197!         sysn  (2,:) = sysn  (3,:) * tms(2,:)
1198!         syysn (2,:) = syysn (3,:) * tms(2,:)
1199!         sxysn (2,:) = sxysn (3,:) * tms(2,:)
1200!         sxa   (2,:) = sxa   (3,:) * tms(2,:)
1201!         sxxa  (2,:) = sxxa  (3,:) * tms(2,:)
1202!         sya   (2,:) = sya   (3,:) * tms(2,:)
1203!         syya  (2,:) = syya  (3,:) * tms(2,:)
1204!         sxya  (2,:) = sxya  (3,:) * tms(2,:)
1205      ENDIF
1206
1207      IF((nbondj == -1).OR.(nbondj == 2)) THEN
1208         frld (:,1) = alpha1 * zai(:,1) + alpha2 * zai(:,2)
1209         hicif(:,1) = alpha1 * zhi(:,1) + alpha2 * zhi(:,2)
1210         hsnif(:,1) = alpha1 * zhs(:,1) + alpha2 * zhs(:,2)
1211! FD add moments for Prather SOM advection
1212         sxice (:,1) = alpha1 * zhim(:,1,1) + alpha2 * zhim(:,2,1)
1213         syice (:,1) = alpha1 * zhim(:,1,2) + alpha2 * zhim(:,2,2)
1214         sxxice(:,1) = alpha1 * zhim(:,1,3) + alpha2 * zhim(:,2,3)
1215         syyice(:,1) = alpha1 * zhim(:,1,4) + alpha2 * zhim(:,2,4)
1216         sxyice(:,1) = alpha1 * zhim(:,1,5) + alpha2 * zhim(:,2,5)
1217         sxa   (:,1) = alpha1 * zaim(:,1,1) + alpha2 * zaim(:,2,1)
1218         sya   (:,1) = alpha1 * zaim(:,1,2) + alpha2 * zaim(:,2,2)
1219         sxxa  (:,1) = alpha1 * zaim(:,1,3) + alpha2 * zaim(:,2,3)
1220         syya  (:,1) = alpha1 * zaim(:,1,4) + alpha2 * zaim(:,2,4)
1221         sxya  (:,1) = alpha1 * zaim(:,1,5) + alpha2 * zaim(:,2,5)
1222
1223         DO ji=1,jpi
1224            IF (vmask(ji,2,1).EQ.0.) THEN
1225               frld (ji,2)=frld (ji,1) * tms(ji,2)
1226               hicif(ji,2)=hicif(ji,1) * tms(ji,2)
1227               hsnif(ji,2)=hsnif(ji,1) * tms(ji,2)
1228! FD add moments for Prather SOM advection
1229               sxice (ji,2) = sxice (ji,1) * tms(ji,2)
1230               syice (ji,2) = syice (ji,1) * tms(ji,2)
1231               sxxice(ji,2) = sxxice(ji,1) * tms(ji,2)
1232               syyice(ji,2) = syyice(ji,1) * tms(ji,2)
1233               sxyice(ji,2) = sxyice(ji,1) * tms(ji,2)
1234               sxa   (ji,2) = sxa   (ji,1) * tms(ji,2)
1235               sya   (ji,2) = sya   (ji,1) * tms(ji,2)
1236               sxxa  (ji,2) = sxxa  (ji,1) * tms(ji,2)
1237               syya  (ji,2) = syya  (ji,1) * tms(ji,2)
1238               sxya  (ji,2) = sxya  (ji,1) * tms(ji,2)
1239
1240            ELSE
1241               frld (ji,2)=(alpha4*frld (ji,1)+alpha3*frld (ji,3))*tms(ji,2)
1242               hicif(ji,2)=(alpha4*hicif(ji,1)+alpha3*hicif(ji,3))*tms(ji,2)
1243               hsnif(ji,2)=(alpha4*hsnif(ji,1)+alpha3*hsnif(ji,3))*tms(ji,2)
1244! FD add moments for Prather SOM advection
1245               sxice (ji,2)=(alpha4*sxice (ji,1)+alpha3*sxice (ji,3))*tms(ji,2)     
1246               syice (ji,2)=(alpha4*syice (ji,1)+alpha3*syice (ji,3))*tms(ji,2)     
1247               sxxice(ji,2)=(alpha4*sxxice(ji,1)+alpha3*sxxice(ji,3))*tms(ji,2) 
1248               syyice(ji,2)=(alpha4*syyice(ji,1)+alpha3*syyice(ji,3))*tms(ji,2) 
1249               sxyice(ji,2)=(alpha4*sxyice(ji,1)+alpha3*sxyice(ji,3))*tms(ji,2) 
1250               sxa   (ji,2)=(alpha4*sxa   (ji,1)+alpha3*sxa   (ji,3))*tms(ji,2)     
1251               sya   (ji,2)=(alpha4*sya   (ji,1)+alpha3*sya   (ji,3))*tms(ji,2)     
1252               sxxa  (ji,2)=(alpha4*sxxa  (ji,1)+alpha3*sxxa  (ji,3))*tms(ji,2) 
1253               syya  (ji,2)=(alpha4*syya  (ji,1)+alpha3*syya  (ji,3))*tms(ji,2) 
1254               sxya  (ji,2)=(alpha4*sxya  (ji,1)+alpha3*sxya  (ji,3))*tms(ji,2) 
1255
1256               IF (zvi_v(ji,2) .LT. 0.) THEN
1257              frld (ji,2)=(alpha6*frld (ji,3)+alpha5*frld (ji,1)+alpha7*frld (ji,4))*tms(ji,2)
1258              hicif(ji,2)=(alpha6*hicif(ji,3)+alpha5*hicif(ji,1)+alpha7*hicif(ji,4))*tms(ji,2)
1259              hsnif(ji,2)=(alpha6*hsnif(ji,3)+alpha5*hsnif(ji,1)+alpha7*hsnif(ji,4))*tms(ji,2)
1260! FD add moments for Prather SOM advection
1261              sxice (ji,2)=(alpha6*sxice (ji,3)+alpha5*sxice (ji,1)+alpha7*sxice (ji,4))*tms(ji,2)
1262              syice (ji,2)=(alpha6*syice (ji,3)+alpha5*syice (ji,1)+alpha7*syice (ji,4))*tms(ji,2)
1263              sxxice(ji,2)=(alpha6*sxxice(ji,3)+alpha5*sxxice(ji,1)+alpha7*sxxice(ji,4))*tms(ji,2)
1264              syyice(ji,2)=(alpha6*syyice(ji,3)+alpha5*syyice(ji,1)+alpha7*syyice(ji,4))*tms(ji,2)
1265              sxyice(ji,2)=(alpha6*sxyice(ji,3)+alpha5*sxyice(ji,1)+alpha7*sxyice(ji,4))*tms(ji,2)
1266              sxa   (ji,2)=(alpha6*sxa   (ji,3)+alpha5*sxa   (ji,1)+alpha7*sxa   (ji,4))*tms(ji,2)
1267              sya   (ji,2)=(alpha6*sya   (ji,3)+alpha5*sya   (ji,1)+alpha7*sya   (ji,4))*tms(ji,2)
1268              sxxa  (ji,2)=(alpha6*sxxa  (ji,3)+alpha5*sxxa  (ji,1)+alpha7*sxxa  (ji,4))*tms(ji,2)
1269              syya  (ji,2)=(alpha6*syya  (ji,3)+alpha5*syya  (ji,1)+alpha7*syya  (ji,4))*tms(ji,2)
1270              sxya  (ji,2)=(alpha6*sxya  (ji,3)+alpha5*sxya  (ji,1)+alpha7*sxya  (ji,4))*tms(ji,2)
1271               ENDIF
1272            ENDIF
1273         END DO
1274!         sxice (:,2) = sxice (:,3) * tms(:,2)
1275!         sxxice(:,2) = sxxice(:,3) * tms(:,2)
1276!         syice (:,2) = syice (:,3) * tms(:,2)
1277!         syyice(:,2) = syyice(:,3) * tms(:,2)
1278!         sxyice(:,2) = sxyice(:,3) * tms(:,2)
1279!         sxsn  (:,2) = sxsn  (:,3) * tms(:,2)
1280!         sxxsn (:,2) = sxxsn (:,3) * tms(:,2)
1281!         sysn  (:,2) = sysn  (:,3) * tms(:,2)
1282!         syysn (:,2) = syysn (:,3) * tms(:,2)
1283!         sxysn (:,2) = sxysn (:,3) * tms(:,2)
1284!         sxa   (:,2) = sxa   (:,3) * tms(:,2)
1285!         sxxa  (:,2) = sxxa  (:,3) * tms(:,2)
1286!         sya   (:,2) = sya   (:,3) * tms(:,2)
1287!         syya  (:,2) = syya  (:,3) * tms(:,2)
1288!         sxya  (:,2) = sxya  (:,3) * tms(:,2)
1289      ENDIF
1290
1291      frld (:,:) = frld (:,:) * tms(:,:)
1292      hicif(:,:) = hicif(:,:) * tms(:,:)
1293      hsnif(:,:) = hsnif(:,:) * tms(:,:)
1294! FD add moments for Prather SOM advection
1295      sxice (:,:) = sxice (:,:) * tms(:,:)
1296      syice (:,:) = syice (:,:) * tms(:,:)
1297      sxxice(:,:) = sxxice(:,:) * tms(:,:)
1298      syyice(:,:) = syyice(:,:) * tms(:,:)
1299      sxyice(:,:) = sxyice(:,:) * tms(:,:)
1300      sxa   (:,:) = sxa   (:,:) * tms(:,:)
1301      sya   (:,:) = sya   (:,:) * tms(:,:)
1302      sxxa  (:,:) = sxxa  (:,:) * tms(:,:)
1303      syya  (:,:) = syya  (:,:) * tms(:,:)
1304      sxya  (:,:) = sxya  (:,:) * tms(:,:)
1305
1306! FD debug
1307if(lwp) write(*,*) 'interp'
1308if(lwp) write(*,'(10(f10.8,1x))') hicif(70,jpj-9:jpj)
1309
1310   END SUBROUTINE Agrif_tra_ice
1311# endif
1312   
1313
1314   SUBROUTINE interpu(tabres,i1,i2,j1,j2,k1,k2)
1315      !!---------------------------------------------
1316      !!   *** ROUTINE interpu ***
1317      !!---------------------------------------------
1318#  include "domzgr_substitute.h90"   
1319   
1320      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1321      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
1322
1323      INTEGER :: ji,jj,jk
1324
1325      DO jk=k1,k2
1326         DO jj=j1,j2
1327            DO ji=i1,i2
1328               tabres(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk)
1329#if ! defined key_zco
1330               tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3u(ji,jj,jk)
1331#endif
1332            END DO
1333         END DO
1334      END DO
1335   END SUBROUTINE interpu
1336
1337# if   defined key_dynspg_ts   ||  defined key_vvl   ||  defined key_esopa
1338   SUBROUTINE interpu2d2(tabres,i1,i2,j1,j2)
1339      !!---------------------------------------------
1340      !!   *** ROUTINE interpu ***
1341      !!---------------------------------------------
1342#  include "domzgr_substitute.h90"   
1343   
1344      INTEGER, INTENT(in) :: i1,i2,j1,j2
1345      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1346
1347      INTEGER :: ji,jj
1348
1349         DO jj=j1,j2
1350            DO ji=i1,i2
1351               tabres(ji,jj) = e2u(ji,jj) * un_b(ji,jj)
1352            END DO
1353         END DO
1354   END SUBROUTINE interpu2d2
1355# endif
1356
1357# if   defined key_ice_lim
1358   SUBROUTINE interpu2di(tabres,i1,i2,j1,j2)
1359      !!---------------------------------------------
1360      !!   *** ROUTINE interpu ***
1361      !!---------------------------------------------
1362#  include "domzgr_substitute.h90"   
1363   
1364      INTEGER, INTENT(in) :: i1,i2,j1,j2
1365      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1366
1367      INTEGER :: ji,jj
1368
1369! FD debug         DO jj=j1+1,j2
1370!            DO ji=i1+1,i2
1371!               tabres(ji,jj) = e2f(ji-1,jj-1) * u_ice(ji,jj)
1372         DO jj=MAX(j1,2),j2
1373            DO ji=MAX(i1,2),i2
1374               tabres(ji,jj) = e2f(ji-1,jj-1) * u_ice(ji,jj)
1375            END DO
1376         END DO
1377   END SUBROUTINE interpu2di
1378# endif
1379
1380   SUBROUTINE interpu2d(tabres,i1,i2,j1,j2)
1381      !!---------------------------------------------
1382      !!   *** ROUTINE interpu2d ***
1383      !!---------------------------------------------
1384
1385      INTEGER, INTENT(in) :: i1,i2,j1,j2
1386      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1387
1388      INTEGER :: ji,jj
1389
1390      DO jj=j1,j2
1391! FD         DO ji=i1,i2 ! out of bound error here
1392         DO ji=i1,i2-1
1393            tabres(ji,jj) = e2u(ji,jj) * ((gcx(ji+1,jj) - gcx(ji,jj))/e1u(ji,jj)) &
1394               * umask(ji,jj,1)
1395         END DO
1396      END DO
1397
1398   END SUBROUTINE interpu2d
1399
1400   SUBROUTINE interpv(tabres,i1,i2,j1,j2,k1,k2)
1401      !!---------------------------------------------
1402      !!   *** ROUTINE interpv ***
1403      !!---------------------------------------------
1404#  include "domzgr_substitute.h90" 
1405     
1406      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1407      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
1408
1409      INTEGER :: ji, jj, jk
1410
1411      DO jk=k1,k2
1412         DO jj=j1,j2
1413            DO ji=i1,i2
1414               tabres(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk)
1415#if ! defined key_zco
1416               tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3v(ji,jj,jk)
1417#endif           
1418            END DO
1419         END DO
1420      END DO
1421
1422   END SUBROUTINE interpv
1423
1424# if   defined key_dynspg_ts   ||  defined key_vvl   ||  defined key_esopa
1425   SUBROUTINE interpv2d2(tabres,i1,i2,j1,j2)
1426      !!---------------------------------------------
1427      !!   *** ROUTINE interpv ***
1428      !!---------------------------------------------
1429#  include "domzgr_substitute.h90" 
1430     
1431      INTEGER, INTENT(in) :: i1,i2,j1,j2
1432      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1433
1434      INTEGER :: ji, jj
1435
1436         DO jj=j1,j2
1437            DO ji=i1,i2
1438               tabres(ji,jj) = e1v(ji,jj) * vn_b(ji,jj)
1439            END DO
1440         END DO
1441
1442   END SUBROUTINE interpv2d2
1443# endif
1444
1445# if   defined key_ice_lim
1446   SUBROUTINE interpv2di(tabres,i1,i2,j1,j2)
1447      !!---------------------------------------------
1448      !!   *** ROUTINE interpv ***
1449      !!---------------------------------------------
1450#  include "domzgr_substitute.h90" 
1451     
1452      INTEGER, INTENT(in) :: i1,i2,j1,j2
1453      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1454
1455      INTEGER :: ji, jj
1456
1457! FD debug         DO jj=j1+1,j2
1458!            DO ji=i1+1,i2
1459!               tabres(ji,jj) = e1f(ji-1,jj-1) * v_ice(ji,jj)
1460         DO jj=MAX(j1,2),j2
1461            DO ji=MAX(i1,2),i2
1462               tabres(ji,jj) = e1f(ji-1,jj-1) * v_ice(ji,jj)
1463            END DO
1464         END DO
1465
1466   END SUBROUTINE interpv2di
1467# endif
1468
1469   SUBROUTINE interpv2d(tabres,i1,i2,j1,j2)
1470      !!---------------------------------------------
1471      !!   *** ROUTINE interpv2d ***
1472      !!---------------------------------------------
1473
1474      INTEGER, INTENT(in) :: i1,i2,j1,j2
1475      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1476
1477      INTEGER :: ji,jj
1478
1479! FD      DO jj=j1,j2 ! out of bound error here
1480      DO jj=j1,j2-1
1481         DO ji=i1,i2
1482            tabres(ji,jj) = e1v(ji,jj) * ((gcx(ji,jj+1) - gcx(ji,jj))/e2v(ji,jj)) &
1483               * vmask(ji,jj,1)
1484         END DO
1485      END DO
1486
1487   END SUBROUTINE interpv2d
1488
1489# if   defined key_ice_lim
1490   SUBROUTINE interp_sxice (tabres,i1,i2,j1,j2)
1491      !!---------------------------------------------
1492      !!   *** ROUTINE interpv ***
1493      !!---------------------------------------------
1494#  include "domzgr_substitute.h90" 
1495     
1496      INTEGER, INTENT(in) :: i1,i2,j1,j2
1497      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1498
1499      INTEGER :: ji, jj
1500
1501         DO jj=j1,j2
1502            DO ji=i1,i2
1503               tabres(ji,jj) = sxice (ji,jj) /area(ji,jj)
1504            END DO
1505         END DO
1506if (lwp) write(*,*) 'interp sxice',sxice(69,92),sxice(69,92)/area(69,92)
1507
1508   END SUBROUTINE interp_sxice
1509
1510   SUBROUTINE interp_syice (tabres,i1,i2,j1,j2)
1511      !!---------------------------------------------
1512      !!   *** ROUTINE interpv ***
1513      !!---------------------------------------------
1514#  include "domzgr_substitute.h90" 
1515     
1516      INTEGER, INTENT(in) :: i1,i2,j1,j2
1517      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1518
1519      INTEGER :: ji, jj
1520
1521         DO jj=j1,j2
1522            DO ji=i1,i2
1523               tabres(ji,jj) = syice (ji,jj) /area(ji,jj)
1524            END DO
1525         END DO
1526
1527   END SUBROUTINE interp_syice
1528
1529   SUBROUTINE interp_sxxice(tabres,i1,i2,j1,j2)
1530      !!---------------------------------------------
1531      !!   *** ROUTINE interpv ***
1532      !!---------------------------------------------
1533#  include "domzgr_substitute.h90" 
1534     
1535      INTEGER, INTENT(in) :: i1,i2,j1,j2
1536      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1537
1538      INTEGER :: ji, jj
1539
1540         DO jj=j1,j2
1541            DO ji=i1,i2
1542               tabres(ji,jj) = sxxice(ji,jj) /area(ji,jj)
1543            END DO
1544         END DO
1545
1546   END SUBROUTINE interp_sxxice
1547
1548   SUBROUTINE interp_syyice(tabres,i1,i2,j1,j2)
1549      !!---------------------------------------------
1550      !!   *** ROUTINE interpv ***
1551      !!---------------------------------------------
1552#  include "domzgr_substitute.h90" 
1553     
1554      INTEGER, INTENT(in) :: i1,i2,j1,j2
1555      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1556
1557      INTEGER :: ji, jj
1558
1559         DO jj=j1,j2
1560            DO ji=i1,i2
1561               tabres(ji,jj) = syyice(ji,jj) /area(ji,jj)
1562            END DO
1563         END DO
1564
1565   END SUBROUTINE interp_syyice
1566
1567   SUBROUTINE interp_sxyice(tabres,i1,i2,j1,j2)
1568      !!---------------------------------------------
1569      !!   *** ROUTINE interpv ***
1570      !!---------------------------------------------
1571#  include "domzgr_substitute.h90" 
1572     
1573      INTEGER, INTENT(in) :: i1,i2,j1,j2
1574      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1575
1576      INTEGER :: ji, jj
1577
1578         DO jj=j1,j2
1579            DO ji=i1,i2
1580               tabres(ji,jj) = sxyice(ji,jj) /area(ji,jj)
1581            END DO
1582         END DO
1583
1584   END SUBROUTINE interp_sxyice
1585
1586   SUBROUTINE interp_sxa (tabres,i1,i2,j1,j2)
1587      !!---------------------------------------------
1588      !!   *** ROUTINE interpv ***
1589      !!---------------------------------------------
1590#  include "domzgr_substitute.h90" 
1591     
1592      INTEGER, INTENT(in) :: i1,i2,j1,j2
1593      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1594
1595      INTEGER :: ji, jj
1596
1597         DO jj=j1,j2
1598            DO ji=i1,i2
1599               tabres(ji,jj) = sxa (ji,jj) /area(ji,jj)
1600            END DO
1601         END DO
1602
1603   END SUBROUTINE interp_sxa
1604
1605   SUBROUTINE interp_sya (tabres,i1,i2,j1,j2)
1606      !!---------------------------------------------
1607      !!   *** ROUTINE interpv ***
1608      !!---------------------------------------------
1609#  include "domzgr_substitute.h90" 
1610     
1611      INTEGER, INTENT(in) :: i1,i2,j1,j2
1612      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1613
1614      INTEGER :: ji, jj
1615
1616         DO jj=j1,j2
1617            DO ji=i1,i2
1618               tabres(ji,jj) = sya (ji,jj) /area(ji,jj)
1619            END DO
1620         END DO
1621
1622   END SUBROUTINE interp_sya
1623
1624   SUBROUTINE interp_sxxa(tabres,i1,i2,j1,j2)
1625      !!---------------------------------------------
1626      !!   *** ROUTINE interpv ***
1627      !!---------------------------------------------
1628#  include "domzgr_substitute.h90" 
1629     
1630      INTEGER, INTENT(in) :: i1,i2,j1,j2
1631      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1632
1633      INTEGER :: ji, jj
1634
1635         DO jj=j1,j2
1636            DO ji=i1,i2
1637               tabres(ji,jj) = sxxa(ji,jj) /area(ji,jj)
1638            END DO
1639         END DO
1640
1641   END SUBROUTINE interp_sxxa
1642
1643   SUBROUTINE interp_syya(tabres,i1,i2,j1,j2)
1644      !!---------------------------------------------
1645      !!   *** ROUTINE interpv ***
1646      !!---------------------------------------------
1647#  include "domzgr_substitute.h90" 
1648     
1649      INTEGER, INTENT(in) :: i1,i2,j1,j2
1650      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1651
1652      INTEGER :: ji, jj
1653
1654         DO jj=j1,j2
1655            DO ji=i1,i2
1656               tabres(ji,jj) = syya(ji,jj) /area(ji,jj)
1657            END DO
1658         END DO
1659
1660   END SUBROUTINE interp_syya
1661
1662   SUBROUTINE interp_sxya(tabres,i1,i2,j1,j2)
1663      !!---------------------------------------------
1664      !!   *** ROUTINE interpv ***
1665      !!---------------------------------------------
1666#  include "domzgr_substitute.h90" 
1667     
1668      INTEGER, INTENT(in) :: i1,i2,j1,j2
1669      REAL(wp),DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1670
1671      INTEGER :: ji, jj
1672
1673         DO jj=j1,j2
1674            DO ji=i1,i2
1675               tabres(ji,jj) = sxya(ji,jj) /area(ji,jj)
1676            END DO
1677         END DO
1678
1679   END SUBROUTINE interp_sxya
1680# endif
1681
1682#else
1683CONTAINS
1684
1685   SUBROUTINE Agrif_OPA_Interp_empty
1686      !!---------------------------------------------
1687      !!   *** ROUTINE agrif_OPA_Interp_empty ***
1688      !!---------------------------------------------
1689      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
1690   END SUBROUTINE Agrif_OPA_Interp_empty
1691#endif
1692END MODULE agrif_opa_interp
1693