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

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

clean online sea ice diagnostics
  • Property svn:keywords set to Id
File size: 29.2 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 limhdf         ! ice horizontal diffusion
21   USE limvar         !
22   USE limadv_prather ! advection scheme (Prather)
23   USE limadv_umx     ! advection scheme (ultimate-macho)
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
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, one can choose between
60      !!     a) an Ultimate-Macho scheme (whose order is defined by nn_limadv_ord) => nn_limadv=0
61      !!     b) and a second order Prather scheme => nn_limadv=-1
62      !!
63      !! ** action :
64      !!---------------------------------------------------------------------
65      INTEGER, INTENT(in) ::   kt   ! number of iteration
66      !
67      INTEGER  ::   ji, jj, jk, jm, jl, jt  ! dummy loop indices
68      INTEGER  ::   initad                  ! number of sub-timestep for the advection
69      REAL(wp) ::   zcfl , zusnit           !   -      -
70      CHARACTER(len=80) :: cltmp
71      !
72      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
73      REAL(wp) ::    zdv, zda
74      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold, zsmvold 
75      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax, zviold, zvsold
76      ! --- diffusion --- !
77      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhdfptab
78      INTEGER , PARAMETER                    ::   ihdf_vars  = 6 ! Number of variables in which we apply horizontal diffusion
79                                                                 !  inside limtrp for each ice category , not counting the
80                                                                 !  variables corresponding to ice_layers
81      ! --- ultimate macho only --- !
82      REAL(wp)                               ::   zdt
83      REAL(wp), POINTER, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box
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_limdiachk )   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      SELECT CASE ( nn_limadv )
160         
161                       !=============================!
162      CASE ( 0 )       !==  Ultimate-MACHO scheme  ==!                   
163                       !=============================!
164     
165         CALL wrk_alloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
166     
167         IF( kt == nit000 .AND. lwp ) THEN
168            WRITE(numout,*)''
169            WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme'
170            WRITE(numout,*)'~~~~~~~~~~~'
171         ENDIF
172         !
173         zdt = rdt_ice / REAL(initad)
174         
175         ! transport
176         zudy(:,:) = u_ice(:,:) * e2u(:,:)
177         zvdx(:,:) = v_ice(:,:) * e1v(:,:)
178         
179         ! define velocity for advection: u*grad(H)
180         DO jj = 2, jpjm1
181            DO ji = fs_2, fs_jpim1
182               IF    ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
183               ELSEIF( u_ice(ji,jj)                  >  0._wp ) THEN   ;   zcu_box(ji,jj) = u_ice(ji-1,jj)
184               ELSE                                                    ;   zcu_box(ji,jj) = u_ice(ji  ,jj)
185               ENDIF
186               
187               IF    ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
188               ELSEIF( v_ice(ji,jj)                  >  0._wp ) THEN   ;   zcv_box(ji,jj) = v_ice(ji,jj-1)
189               ELSE                                                    ;   zcv_box(ji,jj) = v_ice(ji,jj  )
190               ENDIF
191            END DO
192         END DO
193         
194         ! advection
195         DO jt = 1, initad
196            CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) )       ! Open water area
197            DO jl = 1, jpl
198               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) )      ! Ice area
199               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) )      ! Ice  volume
200               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) )    ! Salt content
201               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) )    ! Age content
202               DO jk = 1, nlay_i
203                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) )   ! Ice  heat content
204               END DO
205               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume
206               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content
207            END DO
208         END DO
209         !
210         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
211         DO jl = 2, jpl
212            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
213         END DO
214         !
215         CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
216         
217                       !=============================!
218      CASE ( -1 )      !==     Prather scheme      ==!                   
219                       !=============================!
220
221         CALL wrk_alloc( jpi,jpj,            zarea )
222         CALL wrk_alloc( jpi,jpj,1,          z0opw )
223         CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
224         CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
225         
226         IF( kt == nit000 .AND. lwp ) THEN
227            WRITE(numout,*)''
228            WRITE(numout,*)'lim_adv_xy : Prather advection scheme'
229            WRITE(numout,*)'~~~~~~~~~~~'
230         ENDIF
231         
232         zarea(:,:) = e12t(:,:)
233         
234         !-------------------------
235         ! transported fields                                       
236         !-------------------------
237         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area
238         DO jl = 1, jpl
239            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
240            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
241            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
242            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
243            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
244            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
245            DO jk = 1, nlay_i
246               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
247            END DO
248         END DO
249         
250         
251         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
252            DO jt = 1, initad
253               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
254                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
255               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
256                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
257               DO jl = 1, jpl
258                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
259                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
260                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
261                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
262                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
263                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
264                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
265                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
266                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
267                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
268                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
269                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
270                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
271                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
272                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
273                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
274                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
275                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
276                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
277                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
278                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
279                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
280                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
281                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
282                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
283                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
284                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
285                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
286                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
287                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
288                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
289                  END DO
290               END DO
291            END DO
292         ELSE
293            DO jt = 1, initad
294               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
295                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
296               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
297                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
298               DO jl = 1, jpl
299                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
300                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
301                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
302                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
303                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
304                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
305                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
306                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
307                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
308                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
309                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
310                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
311                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
312                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
313                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
314                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
315                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
316                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
317                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
318                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
319                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
320                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
321                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
322                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
323                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
324                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
325                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
326                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
327                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
328                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
329                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
330                  END DO
331               END DO
332            END DO
333         ENDIF
334         
335         !-------------------------------------------
336         ! Recover the properties from their contents
337         !-------------------------------------------
338         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
339         DO jl = 1, jpl
340            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
341            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
342            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
343            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
344            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
345            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
346            DO jk = 1, nlay_i
347               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
348            END DO
349         END DO
350         
351         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
352         DO jl = 2, jpl
353            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
354         END DO
355         
356         CALL wrk_dealloc( jpi,jpj,            zarea )
357         CALL wrk_dealloc( jpi,jpj,1,          z0opw )
358         CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
359         CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
360
361      END SELECT
362     
363      !------------------------------!
364      ! Diffusion of Ice fields                 
365      !------------------------------!
366      IF( nn_ahi0 /= -1 .AND. nn_limdyn == 2 ) THEN
367         !
368         ! --- Prepare diffusion for variables with categories --- !
369         !     mask eddy diffusivity coefficient at ocean U- and V-points
370         jm=1
371         DO jl = 1, jpl
372            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
373               DO ji = 1 , fs_jpim1
374                  pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,  jl ) ) ) )   &
375                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj)
376                  pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,  jj,  jl ) ) ) )   &
377                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,  jj+1,jl ) ) ) ) * ahiv(ji,jj)
378               END DO
379            END DO
380
381            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1
382            zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1
383            zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1
384            zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1
385            zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
386            zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
387            ! Sample of adding more variables to apply lim_hdf (ihdf_vars must be increased)
388            !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1 
389            !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1
390            DO jk = 1, nlay_i
391              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
392            END DO
393         END DO
394
395         ! --- Prepare diffusion for open water area --- !
396         !     mask eddy diffusivity coefficient at ocean U- and V-points
397         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
398            DO ji = 1 , fs_jpim1
399               pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
400                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
401               pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
402                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
403            END DO
404         END DO
405         !
406         zhdfptab(:,:,jm)= ato_i  (:,:);
407
408         ! --- Apply diffusion --- !
409         CALL lim_hdf( zhdfptab, ihdf_vars )
410
411         ! --- Recover properties --- !
412         jm=1
413         DO jl = 1, jpl
414            a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
415            v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
416            v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
417            smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
418            oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
419            e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1
420            ! Sample of adding more variables to apply lim_hdf
421            !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
422            !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
423            DO jk = 1, nlay_i
424               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1
425            END DO
426         END DO
427         ato_i  (:,:) = zhdfptab(:,:,jm)
428             
429      ENDIF
430
431      ! --- diags ---
432      DO jj = 1, jpj
433         DO ji = 1, jpi
434            diag_trp_ei (ji,jj) = ( SUM( e_i  (ji,jj,1:nlay_i,:) ) -  zeiold(ji,jj) ) * r1_rdtice
435            diag_trp_es (ji,jj) = ( SUM( e_s  (ji,jj,1:nlay_s,:) ) -  zesold(ji,jj) ) * r1_rdtice
436            diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:)          ) - zsmvold(ji,jj) ) * r1_rdtice
437            diag_trp_vi (ji,jj) =   SUM(   v_i(ji,jj,:)            -  zviold(ji,jj,:) ) * r1_rdtice
438            diag_trp_vs (ji,jj) =   SUM(   v_s(ji,jj,:)            -  zvsold(ji,jj,:) ) * r1_rdtice
439         END DO
440      END DO
441     
442      IF( nn_limdyn == 2) THEN
443
444         ! zap small areas
445         CALL lim_var_zapsmall
446           
447         !--- Thickness correction in case too high --- !
448         DO jl = 1, jpl
449            DO jj = 1, jpj
450               DO ji = 1, jpi
451                 
452                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
453                     
454                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
455                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
456                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
457                     
458                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 
459                     
460                     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. &
461                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
462                       
463                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
464                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
465                       
466                        ! small correction due to *rswitch for a_i
467                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
468                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
469                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
470                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
471                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
472                                               
473                     ENDIF
474                     
475                  ENDIF
476               
477               END DO
478            END DO
479         END DO
480         
481         ! Force the upper limit of ht_i to always be < hi_max (99 m).
482         DO jj = 1, jpj
483            DO ji = 1, jpi
484               rswitch         = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
485               ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
486               a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
487            END DO
488         END DO
489
490      ENDIF
491         
492      !------------------------------------------------------------
493      ! Impose a_i < amax if no ridging/rafting or in mono-category
494      !------------------------------------------------------------
495      !
496      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
497      IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2
498         DO jl = 1, jpl
499            DO jj = 1, jpj
500               DO ji = 1, jpi
501                  rswitch       = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) )
502                  zda           = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp )  &
503                     &                    * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 )
504                  a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda
505               END DO
506            END DO
507         END DO
508      ENDIF
509     
510      ! --- agglomerate variables -----------------
511      vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
512      vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
513      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
514     
515      ! --- open water = 1 if at_i=0 --------------------------------
516      WHERE( at_i == 0._wp ) ato_i = 1._wp 
517     
518      ! conservation test
519      IF( ln_limdiachk )   CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
520       
521      ! -------------------------------------------------
522      ! control prints
523      ! -------------------------------------------------
524      IF( ln_limctl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
525      !
526      CALL wrk_dealloc( jpi,jpj,                            zatold, zeiold, zesold, zsmvold )
527      CALL wrk_dealloc( jpi,jpj,jpl,                        zhimax, zviold, zvsold )
528      CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1, zhdfptab)
529      !
530      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
531      !
532   END SUBROUTINE lim_trp
533
534#else
535   !!----------------------------------------------------------------------
536   !!   Default option         Empty Module                No sea-ice model
537   !!----------------------------------------------------------------------
538CONTAINS
539   SUBROUTINE lim_trp        ! Empty routine
540   END SUBROUTINE lim_trp
541#endif
542
543   !!======================================================================
544END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.