source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limadv_prather.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 4 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

File size: 22.7 KB
Line 
1MODULE limadv_prather 
2   !!======================================================================
3   !!                       ***  MODULE limadv_prather   ***
4   !! LIM sea-ice model : sea-ice advection => Prather scheme
5   !!======================================================================
6   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code
7   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!--------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                     LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_adv_x  : advection of sea ice on x axis
15   !!   lim_adv_y  : advection of sea ice on y axis
16   !!----------------------------------------------------------------------
17   USE dom_oce          ! ocean domain
18   USE ice              ! LIM-3 variables
19   USE lbclnk           ! lateral boundary condition - MPP exchanges
20   USE in_out_manager   ! I/O manager
21   USE prtctl           ! Print control
22   USE lib_mpp          ! MPP library
23   USE wrk_nemo         ! work arrays
24   USE lib_fortran      ! to use key_nosignedzero
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   lim_adv_x   ! called by lim_trp
30   PUBLIC   lim_adv_y   ! called by lim_trp
31
32   !! * Substitutions
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
36   !! $Id: limadv.F90 6746 2016-06-27 17:20:57Z clem $
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE lim_adv_x( pdf, put , pcrh, psm , ps0 ,   &
42      &                  psx, psxx, psy , psyy, psxy )
43      !!---------------------------------------------------------------------
44      !!                **  routine lim_adv_x  **
45      !! 
46      !! ** purpose :   Computes and adds the advection trend to sea-ice
47      !!              variable on i-axis
48      !!
49      !! ** method  :   Uses Prather second order scheme that advects tracers
50      !!              but also theirquadratic forms. The method preserves
51      !!              tracer structures by conserving second order moments.
52      !!
53      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
54      !!--------------------------------------------------------------------
55      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
56      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
57      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
58      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
59      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
60      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
61      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
62      !!
63      INTEGER  ::   ji, jj                               ! dummy loop indices
64      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
65      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
66      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
67      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
68      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
69      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q         !  -      -
70      !---------------------------------------------------------------------
71
72      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
73      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
74
75      ! Limitation of moments.                                           
76
77      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
78
79      DO jj = 1, jpj
80         DO ji = 1, jpi
81            zslpmax = MAX( 0._wp, ps0(ji,jj) )
82            zs1max  = 1.5 * zslpmax
83            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
84            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
85               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
86            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
87
88            ps0 (ji,jj) = zslpmax 
89            psx (ji,jj) = zs1new      * rswitch
90            psxx(ji,jj) = zs2new      * rswitch
91            psy (ji,jj) = psy (ji,jj) * rswitch
92            psyy(ji,jj) = psyy(ji,jj) * rswitch
93            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
94         END DO
95      END DO
96
97      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
98      psm (:,:)  = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
99
100      !  Calculate fluxes and moments between boxes i<-->i+1             
101      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0
102         DO ji = 1, jpi
103            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
104            zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
105            zalfq        =  zalf * zalf
106            zalf1        =  1.0 - zalf
107            zalf1q       =  zalf1 * zalf1
108            !
109            zfm (ji,jj)  =  zalf  *   psm (ji,jj)
110            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  )
111            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
112            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq
113            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )
114            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj)
115            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj)
116
117            !  Readjust moments remaining in the box.
118            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
119            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
120            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
121            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
122            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
123            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
124            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
125         END DO
126      END DO
127
128      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
129         DO ji = 1, fs_jpim1
130            zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
131            zalg  (ji,jj) = zalf
132            zalfq         = zalf * zalf
133            zalf1         = 1.0 - zalf
134            zalg1 (ji,jj) = zalf1
135            zalf1q        = zalf1 * zalf1
136            zalg1q(ji,jj) = zalf1q
137            !
138            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj)
139            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
140            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
141            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq
142            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )
143            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj)
144            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj)
145         END DO
146      END DO
147
148      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box.
149         DO ji = fs_2, fs_jpim1
150            zbt  =       zbet(ji-1,jj)
151            zbt1 = 1.0 - zbet(ji-1,jj)
152            !
153            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
154            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
155            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
156            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
157            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
158            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
159            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
160         END DO
161      END DO
162
163      !   Put the temporary moments into appropriate neighboring boxes.   
164      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
165         DO ji = fs_2, fs_jpim1
166            zbt  =       zbet(ji-1,jj)
167            zbt1 = 1.0 - zbet(ji-1,jj)
168            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
169            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
170            zalf1       = 1.0 - zalf
171            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
172            !
173            ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)
174            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
175            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
176               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
177               &                                                + zbt1 * psxx(ji,jj)
178            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
179               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
180               &                                                + zbt1 * psxy(ji,jj)
181            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
182            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
183         END DO
184      END DO
185
186      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
187         DO ji = fs_2, fs_jpim1
188            zbt  =       zbet(ji,jj)
189            zbt1 = 1.0 - zbet(ji,jj)
190            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
191            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
192            zalf1       = 1.0 - zalf
193            ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
194            !
195            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
196            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
197            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  &
198               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      &
199               &                                      + ( zalf1 - zalf ) * ztemp ) )
200            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  &
201               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
202            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
203            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
204         END DO
205      END DO
206
207      !-- Lateral boundary conditions
208      CALL lbc_lnk_multi( psm , 'T',  1., ps0 , 'T',  1.   &
209         &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
210         &              , psxx, 'T',  1., psyy, 'T',  1.   &
211         &              , psxy, 'T',  1. )
212
213      IF(ln_ctl) THEN
214         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
215         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
216         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
217         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
218      ENDIF
219      !
220      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
221      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
222      !
223   END SUBROUTINE lim_adv_x
224
225
226   SUBROUTINE lim_adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
227      &                  psx, psxx, psy , psyy, psxy )
228      !!---------------------------------------------------------------------
229      !!                **  routine lim_adv_y  **
230      !!           
231      !! ** purpose :   Computes and adds the advection trend to sea-ice
232      !!              variable on y axis
233      !!
234      !! ** method  :   Uses Prather second order scheme that advects tracers
235      !!              but also their quadratic forms. The method preserves
236      !!              tracer structures by conserving second order moments.
237      !!
238      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
239      !!---------------------------------------------------------------------
240      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
241      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
242      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
243      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
244      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
245      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments
246      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
247      !!
248      INTEGER  ::   ji, jj                               ! dummy loop indices
249      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
250      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
251      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
252      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0, zfx , zfy , zbet   ! 2D workspace
253      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm, zfxx, zfyy, zfxy   !  -      -
254      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q     !  -      -
255      !---------------------------------------------------------------------
256
257      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
258      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
259
260      ! Limitation of moments.
261
262      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
263
264      DO jj = 1, jpj
265         DO ji = 1, jpi
266            zslpmax = MAX( 0._wp, ps0(ji,jj) )
267            zs1max  = 1.5 * zslpmax
268            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
269            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
270               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
271            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
272            !
273            ps0 (ji,jj) = zslpmax 
274            psx (ji,jj) = psx (ji,jj) * rswitch
275            psxx(ji,jj) = psxx(ji,jj) * rswitch
276            psy (ji,jj) = zs1new * rswitch
277            psyy(ji,jj) = zs2new * rswitch
278            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
279         END DO
280      END DO
281
282      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
283      psm(:,:)  = MAX(  pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
284
285      !  Calculate fluxes and moments between boxes j<-->j+1             
286      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
287         DO ji = 1, jpi
288            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
289            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
290            zalfq        =  zalf * zalf
291            zalf1        =  1.0 - zalf
292            zalf1q       =  zalf1 * zalf1
293            !
294            zfm (ji,jj)  =  zalf  * psm(ji,jj)
295            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
296            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
297            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
298            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
299            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
300            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
301            !
302            !  Readjust moments remaining in the box.
303            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
304            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
305            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
306            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
307            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
308            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
309            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
310         END DO
311      END DO
312      !
313      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
314         DO ji = 1, jpi
315            zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
316            zalg  (ji,jj) = zalf
317            zalfq         = zalf * zalf
318            zalf1         = 1.0 - zalf
319            zalg1 (ji,jj) = zalf1
320            zalf1q        = zalf1 * zalf1
321            zalg1q(ji,jj) = zalf1q
322            !
323            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
324            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
325            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
326            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
327            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
328            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
329            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
330         END DO
331      END DO
332
333      !  Readjust moments remaining in the box.
334      DO jj = 2, jpj
335         DO ji = 1, jpi
336            zbt  =         zbet(ji,jj-1)
337            zbt1 = ( 1.0 - zbet(ji,jj-1) )
338            !
339            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
340            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
341            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
342            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
343            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
344            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
345            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
346         END DO
347      END DO
348
349      !   Put the temporary moments into appropriate neighboring boxes.   
350      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
351         DO ji = 1, jpi
352            zbt  =         zbet(ji,jj-1)
353            zbt1 = ( 1.0 - zbet(ji,jj-1) )
354            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
355            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
356            zalf1       = 1.0 - zalf
357            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
358            !
359            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
360            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
361               &                                               + zbt1 * psy(ji,jj) 
362            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
363               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
364               &                                               + zbt1 * psyy(ji,jj)
365            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
366               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
367               &                                                + zbt1 * psxy(ji,jj)
368            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
369            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
370         END DO
371      END DO
372
373      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
374         DO ji = 1, jpi
375            zbt  =         zbet(ji,jj)
376            zbt1 = ( 1.0 - zbet(ji,jj) )
377            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
378            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
379            zalf1       = 1.0 - zalf
380            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
381            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
382            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
383            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
384               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
385               &                                         + ( zalf1 - zalf ) * ztemp )                                )
386            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
387               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
388            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
389            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
390         END DO
391      END DO
392
393      !-- Lateral boundary conditions
394      CALL lbc_lnk_multi( psm , 'T',  1.,  ps0 , 'T',  1.   &
395         &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes
396         &              , psxx, 'T',  1.,  psyy, 'T',  1.   &
397         &              , psxy, 'T',  1. )
398
399      IF(ln_ctl) THEN
400         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
401         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
402         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
403         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
404      ENDIF
405      !
406      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
407      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
408      !
409   END SUBROUTINE lim_adv_y
410
411#else
412   !!----------------------------------------------------------------------
413   !!   Default option            Dummy module         NO LIM sea-ice model
414   !!----------------------------------------------------------------------
415#endif
416
417   !!======================================================================
418END MODULE limadv_prather
Note: See TracBrowser for help on using the repository browser.