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.
limtrp.F90 in branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 6584

Last change on this file since 6584 was 6584, checked in by clem, 8 years ago

LIM3 and Agrif compatibility
  • Property svn:keywords set to Id
File size: 29.0 KB
Line 
1MODULE limtrp
2   !!======================================================================
3   !!                       ***  MODULE limtrp   ***
4   !! LIM transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6   !! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code
7   !!            3.0  ! 2005-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
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_trp      : advection/diffusion process of sea ice
15   !!----------------------------------------------------------------------
16   USE phycst         ! physical constant
17   USE dom_oce        ! ocean domain
18   USE sbc_oce        ! ocean surface boundary condition
19   USE dom_ice        ! ice domain
20   USE ice            ! ice variables
21   USE limadv         ! ice advection
22   USE limhdf         ! ice horizontal diffusion
23   USE limvar         !
24   !
25   USE in_out_manager ! I/O manager
26   USE lbclnk         ! lateral boundary conditions -- MPP exchanges
27   USE lib_mpp        ! MPP library
28   USE wrk_nemo       ! work arrays
29   USE prtctl         ! Print control
30   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
31   USE timing         ! Timing
32   USE limcons        ! conservation tests
33   USE limctl         ! control prints
34   USE limadv_umx     ! advection scheme
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   lim_trp    ! called by sbcice_lim
40
41   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2 
42
43   !! * Substitution
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE lim_trp( kt ) 
53      !!-------------------------------------------------------------------
54      !!                   ***  ROUTINE lim_trp ***
55      !!                   
56      !! ** purpose : advection/diffusion process of sea ice
57      !!
58      !! ** method  : variables included in the process are scalar,   
59      !!     other values are considered as second order.
60      !!     For advection, a second order Prather scheme is used. 
61      !!
62      !! ** action :
63      !!---------------------------------------------------------------------
64      INTEGER, INTENT(in) ::   kt   ! number of iteration
65      !
66      INTEGER  ::   ji, jj, jk, jm, jl, jt  ! dummy loop indices
67      INTEGER  ::   initad                  ! number of sub-timestep for the advection
68      REAL(wp) ::   zcfl , zusnit           !   -      -
69      CHARACTER(len=80) :: cltmp
70      !
71      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
72      REAL(wp) ::    zdv, zda
73      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold, zsmvold 
74      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax, zviold, zvsold
75      ! --- diffusion --- !
76      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhdfptab
77      INTEGER , PARAMETER                    ::   ihdf_vars  = 6 ! Number of variables in which we apply horizontal diffusion
78                                                                 !  inside limtrp for each ice category , not counting the
79                                                                 !  variables corresponding to ice_layers
80      ! --- ultimate macho only --- !
81      REAL(wp)                               ::   zdt
82      LOGICAL                                ::   lcon
83      REAL(wp), POINTER, DIMENSION(:,:)      ::   ze, zu_trp, zv_trp, z1_v, zudy, zvdx
84      ! --- prather only --- !
85      REAL(wp), POINTER, DIMENSION(:,:)      ::   zarea
86      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
87      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
88      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
89      !!
90      !!---------------------------------------------------------------------
91      IF( nn_timing == 1 )  CALL timing_start('limtrp')
92
93      CALL wrk_alloc( jpi,jpj,                            zatold, zeiold, zesold, zsmvold )
94      CALL wrk_alloc( jpi,jpj,jpl,                        zhimax, zviold, zvsold )
95      CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1, zhdfptab)
96 
97      IF( kt == nit000 .AND. lwp ) THEN
98         WRITE(numout,*)''
99         WRITE(numout,*)'limtrp'
100         WRITE(numout,*)'~~~~~~'
101         ncfl = 0                ! nb of time step with CFL > 1/2
102      ENDIF
103     
104      CALL lim_var_agg( 1 ) ! integrated values + ato_i
105
106      !-------------------------------------!
107      !   Advection of sea ice properties   !
108      !-------------------------------------!
109
110      ! conservation test
111      IF( ln_limdiahsb )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
112     
113      ! store old values for diag
114      zviold = v_i
115      zvsold = v_s
116      zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 )
117      zeiold (:,:) = et_i
118      zesold (:,:) = et_s 
119
120      !--- Thickness correction init. --- !
121      zatold(:,:) = at_i
122      DO jl = 1, jpl
123         DO jj = 1, jpj
124            DO ji = 1, jpi
125               rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
126               ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
127               ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
128            END DO
129         END DO
130      END DO
131      ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- !
132      zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
133      DO jl = 1, jpl
134         DO jj = 2, jpjm1
135            DO ji = 2, jpim1
136               zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )
137            END DO
138         END DO
139         CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
140      END DO
141         
142      ! --- If ice drift field is too fast, use an appropriate time step for advection --- !       
143      zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )    ! CFL test for stability
144      zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
145      IF(lk_mpp )   CALL mpp_max( zcfl )
146     
147      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
148      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
149      ENDIF
150     
151!!      IF( zcfl > 0.5_wp .AND. lwp ) THEN
152!!         ncfl = ncfl + 1
153!!         IF( ncfl > 0 ) THEN   
154!!            WRITE(cltmp,'(i6.1)') ncfl
155!!            CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
156!!         ENDIF
157!!      ENDIF
158
159#if defined key_limumx
160      !=============================!
161      !==  Ultimate-MACHO scheme  ==!
162      !=============================!
163     
164      CALL wrk_alloc( jpi,jpj, ze, zu_trp, zv_trp, z1_v, zudy, zvdx )
165     
166      IF( kt == nit000 .AND. lwp ) THEN
167         WRITE(numout,*)''
168         WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme'
169         WRITE(numout,*)'~~~~~~~~~~~'
170      ENDIF
171      !
172      zdt = rdt_ice / REAL(initad)
173     
174      ! transport
175      zudy(:,:) = u_ice(:,:) * e2u(:,:)
176      zvdx(:,:) = v_ice(:,:) * e1v(:,:)
177      !
178      DO jt = 1, initad
179         lcon = .TRUE. 
180!!!            lcon = .false.
181         CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, ato_i(:,:), ato_i(:,:) )                                       ! Open water area
182         !
183         DO jl = 1, jpl
184            WHERE( v_i(:,:,jl) /= 0._wp )   ;   z1_v(:,:) = 1._wp / v_i(:,:,jl)
185            ELSEWHERE                       ;   z1_v(:,:) = 0._wp
186            END WHERE
187            !
188            lcon = .TRUE. 
189!!!               lcon = .false.
190            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, a_i(:,:,jl), a_i(:,:,jl) )                      ! Ice area
191            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, v_i(:,:,jl), v_i(:,:,jl), zu_trp, zv_trp )      ! Ice  volume
192            !
193            lcon = .FALSE. 
194            ze(:,:) = smv_i(:,:,jl) * z1_v(:,:)
195            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, smv_i(:,:,jl) )                           ! Salt content
196            !
197!!!check that               ze(:,:) = oa_i (:,:,jl) * z1_v(:,:)
198!!!check that               CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, oa_i (:,:,jl) )                           ! Age content
199            !
200            zu_trp(:,:) = zu_trp(:,:) * r1_nlay_i
201            zv_trp(:,:) = zv_trp(:,:) * r1_nlay_i
202            z1_v  (:,:) = z1_v  (:,:) * REAL( nlay_i, wp )
203            DO jk = 1, nlay_i
204               ze (:,:) = e_i(:,:,jk,jl) * z1_v(:,:)
205               CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, e_i(:,:,jk,jl) )                                  ! Ice  heat content
206            END DO
207            !
208            WHERE( v_s(:,:,jl) /= 0._wp )   ;   z1_v(:,:) = 1._wp / v_s(:,:,jl)
209            ELSEWHERE                       ;   z1_v(:,:) = 0._wp
210            END WHERE
211            !
212            lcon = .TRUE. 
213!!!               lcon = .false.
214            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, v_s(:,:,jl), v_s(:,:,jl), zu_trp, zv_trp )      ! Snow volume
215            !
216            lcon = .FALSE. 
217            ze (:,:) = e_s(:,:,1,jl) * z1_v(:,:)
218            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, e_s(:,:,1,jl) )                                     ! Snow heat content
219            !
220         END DO
221      END DO
222      !
223      at_i(:,:) = a_i(:,:,1)      ! total ice fraction
224      DO jl = 2, jpl
225         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
226      END DO
227      !
228      CALL wrk_dealloc( jpi,jpj, ze, zu_trp, zv_trp, z1_v, zudy, zvdx )
229       
230#else
231      !=============================!
232      !==      Prather scheme     ==!
233      !=============================!
234
235      CALL wrk_alloc( jpi,jpj,            zarea )
236      CALL wrk_alloc( jpi,jpj,1,          z0opw )
237      CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
238      CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
239     
240      IF( kt == nit000 .AND. lwp ) THEN
241         WRITE(numout,*)''
242         WRITE(numout,*)'lim_adv_xy : Prather advection scheme'
243         WRITE(numout,*)'~~~~~~~~~~~'
244      ENDIF
245
246      zarea(:,:) = e12t(:,:)
247     
248      !-------------------------
249      ! transported fields                                       
250      !-------------------------
251      z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area
252      DO jl = 1, jpl
253         z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
254         z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
255         z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
256         z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
257         z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
258         z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
259         DO jk = 1, nlay_i
260            z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
261         END DO
262      END DO
263
264     
265      IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
266         DO jt = 1, initad
267            CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
268               &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
269            CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
270               &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
271            DO jl = 1, jpl
272               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
273                  &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
274               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
275                  &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
276               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
277                  &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
278               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
279                  &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
280               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
281                  &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
282               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
283                  &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
284               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
285                  &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
286               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
287                  &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
288               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
289                  &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
290               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
291                  &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
292               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
293                  &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
294               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
295                  &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
296               DO jk = 1, nlay_i                                                                !--- ice heat contents ---
297                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
298                     &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
299                     &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
300                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
301                     &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
302                     &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
303               END DO
304            END DO
305         END DO
306      ELSE
307         DO jt = 1, initad
308            CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
309               &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
310            CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
311               &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
312            DO jl = 1, jpl
313               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
314                  &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
315               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
316                  &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
317               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
318                  &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
319               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
320                  &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
321               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
322                  &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
323               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
324                  &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
325               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
326                  &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
327               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
328                  &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
329               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
330                  &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
331               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
332                  &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
333               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
334                  &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
335               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
336                  &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
337               DO jk = 1, nlay_i                                                           !--- ice heat contents ---
338                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
339                     &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
340                     &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
341                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
342                     &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
343                     &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
344               END DO
345            END DO
346         END DO
347      ENDIF
348     
349      !-------------------------------------------
350      ! Recover the properties from their contents
351      !-------------------------------------------
352      ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
353      DO jl = 1, jpl
354         v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
355         v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
356         smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
357         oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
358         a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
359         e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
360         DO jk = 1, nlay_i
361            e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
362         END DO
363      END DO
364     
365      at_i(:,:) = a_i(:,:,1)      ! total ice fraction
366      DO jl = 2, jpl
367         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
368      END DO
369
370      CALL wrk_dealloc( jpi,jpj,            zarea )
371      CALL wrk_dealloc( jpi,jpj,1,          z0opw )
372      CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
373      CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
374     
375#endif
376
377      !------------------------------!
378      ! Diffusion of Ice fields                 
379      !------------------------------!
380      IF( nn_ahi0 /= -1 .AND. nn_limdyn == 2 ) THEN
381         !
382         ! --- Prepare diffusion for variables with categories --- !
383         !     mask eddy diffusivity coefficient at ocean U- and V-points
384         jm=1
385         DO jl = 1, jpl
386            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
387               DO ji = 1 , fs_jpim1
388                  pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,  jl ) ) ) )   &
389                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj)
390                  pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,  jj,  jl ) ) ) )   &
391                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,  jj+1,jl ) ) ) ) * ahiv(ji,jj)
392               END DO
393            END DO
394
395            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1
396            zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1
397            zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1
398            zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1
399            zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
400            zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
401            ! Sample of adding more variables to apply lim_hdf (ihdf_vars must be increased)
402            !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1 
403            !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1
404            DO jk = 1, nlay_i
405              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
406            END DO
407         END DO
408
409         ! --- Prepare diffusion for open water area --- !
410         !     mask eddy diffusivity coefficient at ocean U- and V-points
411         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
412            DO ji = 1 , fs_jpim1
413               pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
414                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
415               pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
416                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
417            END DO
418         END DO
419         !
420         zhdfptab(:,:,jm)= ato_i  (:,:);
421
422         ! --- Apply diffusion --- !
423         CALL lim_hdf( zhdfptab, ihdf_vars )
424
425         ! --- Recover properties --- !
426         jm=1
427         DO jl = 1, jpl
428            a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
429            v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
430            v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
431            smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
432            oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
433            e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1
434            ! Sample of adding more variables to apply lim_hdf
435            !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
436            !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
437            DO jk = 1, nlay_i
438               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1
439            END DO
440         END DO
441         ato_i  (:,:) = zhdfptab(:,:,jm)
442             
443      ENDIF
444
445      ! --- diags ---
446      DO jj = 1, jpj
447         DO ji = 1, jpi
448            diag_trp_ei (ji,jj) = ( SUM( e_i  (ji,jj,1:nlay_i,:) ) -  zeiold(ji,jj) ) * r1_rdtice
449            diag_trp_es (ji,jj) = ( SUM( e_s  (ji,jj,1:nlay_s,:) ) -  zesold(ji,jj) ) * r1_rdtice
450            diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:)          ) - zsmvold(ji,jj) ) * r1_rdtice
451            diag_trp_vi (ji,jj) =   SUM(   v_i(ji,jj,:)            -  zviold(ji,jj,:) ) * r1_rdtice
452            diag_trp_vs (ji,jj) =   SUM(   v_s(ji,jj,:)            -  zvsold(ji,jj,:) ) * r1_rdtice
453         END DO
454      END DO
455     
456      IF( nn_limdyn == 2) THEN
457
458         ! zap small areas
459         CALL lim_var_zapsmall
460           
461         !--- Thickness correction in case too high --- !
462         DO jl = 1, jpl
463            DO jj = 1, jpj
464               DO ji = 1, jpi
465                 
466                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
467                     
468                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
469                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
470                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
471                     
472                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 
473                     
474                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
475                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
476                       
477                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
478                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
479                       
480                        ! small correction due to *rswitch for a_i
481                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
482                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
483                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
484                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
485                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
486                                               
487                     ENDIF
488                     
489                  ENDIF
490               
491               END DO
492            END DO
493         END DO
494         
495         ! Force the upper limit of ht_i to always be < hi_max (99 m).
496         DO jj = 1, jpj
497            DO ji = 1, jpi
498               rswitch         = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
499               ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
500               a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
501            END DO
502         END DO
503
504      ENDIF
505         
506      !------------------------------------------------------------
507      ! Impose a_i < amax if no ridging/rafting or in mono-category
508      !------------------------------------------------------------
509      !
510      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
511      IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2
512         DO jl = 1, jpl
513            DO jj = 1, jpj
514               DO ji = 1, jpi
515                  rswitch       = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) )
516                  zda           = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp )  &
517                     &                    * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 )
518                  a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda
519               END DO
520            END DO
521         END DO
522      ENDIF
523     
524      ! --- agglomerate variables -----------------
525      vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
526      vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
527      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
528     
529      ! --- open water = 1 if at_i=0 --------------------------------
530      WHERE( at_i == 0._wp ) ato_i = 1._wp 
531     
532      ! conservation test
533      IF( ln_limdiahsb )   CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
534       
535      ! -------------------------------------------------
536      ! control prints
537      ! -------------------------------------------------
538      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
539      !
540      CALL wrk_dealloc( jpi,jpj,                            zatold, zeiold, zesold, zsmvold )
541      CALL wrk_dealloc( jpi,jpj,jpl,                        zhimax, zviold, zvsold )
542      CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1, zhdfptab)
543      !
544      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
545      !
546   END SUBROUTINE lim_trp
547
548#else
549   !!----------------------------------------------------------------------
550   !!   Default option         Empty Module                No sea-ice model
551   !!----------------------------------------------------------------------
552CONTAINS
553   SUBROUTINE lim_trp        ! Empty routine
554   END SUBROUTINE lim_trp
555#endif
556
557   !!======================================================================
558END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.