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 @ 6746

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

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