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/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 9270

Last change on this file since 9270 was 7993, checked in by frrh, 7 years ago

Merge in missing revisions 6428:2477 inclusive and 6482 from nemo_v3_6_STABLE
branch. In ptic, this includes the fix for restartability of runoff fields in coupled
models. Evolution of coupled models will therefor be affected.

These changes donot affect evolution of the current stand-alone NEMO-CICE GO6
standard configuration.

Work and testing documented in Met Office GMED ticket 320.

File size: 29.4 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
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), POINTER, DIMENSION(:,:)      ::   zsm
71      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
72      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
73      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
74      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold, zsmvold  ! old ice volume...
75      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness
76      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies
77      REAL(wp), POINTER, DIMENSION(:,:,:)             ::   zhdfptab
78      REAL(wp) ::    zdv, zvi, zvs, zsmv, zes, zei
79      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
80      !!---------------------------------------------------------------------
81      INTEGER                                ::  ihdf_vars  = 6  !!Number of variables in which we apply horizontal diffusion
82                                                                   !!  inside limtrp for each ice category , not counting the
83                                                                   !!  variables corresponding to ice_layers
84      !!---------------------------------------------------------------------
85      IF( nn_timing == 1 )  CALL timing_start('limtrp')
86
87      CALL wrk_alloc( jpi,jpj,            zsm, zatold, zeiold, zesold )
88      CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
89      CALL wrk_alloc( jpi,jpj,1,          z0opw )
90      CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
91      CALL wrk_alloc( jpi,jpj,jpl,        zhimax, zviold, zvsold, zsmvold )
92      CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1,zhdfptab)
93
94      IF( numit == nstart .AND. lwp ) THEN
95         WRITE(numout,*)
96         IF( ln_limdyn ) THEN   ;   WRITE(numout,*) 'lim_trp : Ice transport '
97         ELSE                   ;   WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
98         ENDIF
99         WRITE(numout,*) '~~~~~~~~~~~~'
100         ncfl = 0                ! nb of time step with CFL > 1/2
101      ENDIF
102
103      zsm(:,:) = e12t(:,:)
104     
105      !                             !-------------------------------------!
106      IF( ln_limdyn ) THEN          !   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         ! mass and salt flux init
113         zviold(:,:,:)  = v_i(:,:,:)
114         zvsold(:,:,:)  = v_s(:,:,:)
115         zsmvold(:,:,:) = smv_i(:,:,:)
116         zeiold(:,:)    = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 
117         zesold(:,:)    = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 
118
119         !--- Thickness correction init. -------------------------------
120         zatold(:,:) = SUM( a_i(:,:,:), dim=3 )
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         !---------------------------------------------------------------------
131         ! Record max of the surrounding ice thicknesses for correction
132         ! in case advection creates ice too thick.
133         !---------------------------------------------------------------------
134         zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
135         DO jl = 1, jpl
136            DO jj = 2, jpjm1
137               DO ji = 2, jpim1
138                  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) )
139               END DO
140            END DO
141            CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
142         END DO
143         
144         !=============================!
145         !==      Prather scheme     ==!
146         !=============================!
147
148         ! If ice drift field is too fast, use an appropriate time step for advection.         
149         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )         ! CFL test for stability
150         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
151         IF(lk_mpp )   CALL mpp_max( zcfl )
152
153         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
154         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
155         ENDIF
156
157         IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1
158!!         IF( lwp ) THEN
159!!            IF( ncfl > 0 ) THEN   
160!!               WRITE(cltmp,'(i6.1)') ncfl
161!!               CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
162!!            ELSE
163!!            !  WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 '
164!!            ENDIF
165!!         ENDIF
166
167         !-------------------------
168         ! transported fields                                       
169         !-------------------------
170         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area
171         DO jl = 1, jpl
172            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
173            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
174            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
175            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
176            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
177            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
178           DO jk = 1, nlay_i
179               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
180            END DO
181         END DO
182
183
184         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
185            DO jt = 1, initad
186               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
187                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
188               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
189                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
190               DO jl = 1, jpl
191                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
192                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
193                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
194                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
195                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
196                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
197                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
198                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
199                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
200                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
201                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
202                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
203                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
204                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
205                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
206                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
207                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
208                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
209                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
210                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
211                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
212                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
213                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
214                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
215                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
216                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
217                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
218                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
219                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
220                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
221                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
222                  END DO
223               END DO
224            END DO
225         ELSE
226            DO jt = 1, initad
227               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
228                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
229               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
230                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
231               DO jl = 1, jpl
232                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
233                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
234                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
235                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
236                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
237                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
238                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
239                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
240                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
241                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
242                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
243                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
244                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
245                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
246                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
247                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
248                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
249                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
250                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &
251                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
252                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
253                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
254                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
255                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
256                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
257                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
258                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
259                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
260                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
261                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
262                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
263                  END DO
264               END DO
265            END DO
266         ENDIF
267
268         !-------------------------------------------
269         ! Recover the properties from their contents
270         !-------------------------------------------
271         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
272         DO jl = 1, jpl
273            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
274            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
275            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
276            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
277            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
278            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
279            DO jk = 1, nlay_i
280               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
281            END DO
282         END DO
283
284         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
285         DO jl = 2, jpl
286            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
287         END DO
288
289         !------------------------------------------------------------------------------!
290         ! Diffusion of Ice fields                 
291         !------------------------------------------------------------------------------!
292         !------------------------------------
293         !  Diffusion of other ice variables
294         !------------------------------------
295         jm=1
296         DO jl = 1, jpl
297         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
298         !   DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
299         !      DO ji = 1 , fs_jpim1   ! vector opt.
300         !         pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   &
301         !            &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
302         !         pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   &
303         !            &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
304         !      END DO
305         !   END DO
306            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
307               DO ji = 1 , fs_jpim1   ! vector opt.
308                  pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,  jl ) ) ) )   &
309                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj)
310                  pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,  jj,  jl ) ) ) )   &
311                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,  jj+1,jl ) ) ) ) * ahiv(ji,jj)
312               END DO
313            END DO
314
315            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1   
316            zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1
317            zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1 
318            zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1
319            zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
320            zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
321         ! Sample of adding more variables to apply lim_hdf using lim_hdf optimization---
322         !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1 
323         !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1
324         !
325         ! and in this example the parameter ihdf_vars musb be changed to 8 (necessary for allocation)
326         !----------------------------------------------------------------------------------------
327            DO jk = 1, nlay_i
328              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
329            END DO
330         END DO
331         !
332         !--------------------------------
333         !  diffusion of open water area
334         !--------------------------------
335         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
336         !DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
337         !   DO ji = 1 , fs_jpim1   ! vector opt.
338         !      pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
339         !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
340         !      pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
341         !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
342         !   END DO
343         !END DO
344         
345         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
346            DO ji = 1 , fs_jpim1   ! vector opt.
347               pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
348                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
349               pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
350                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
351            END DO
352         END DO
353         !
354         zhdfptab(:,:,jm)= ato_i  (:,:);
355         CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i) 
356
357         jm=1
358         DO jl = 1, jpl
359            a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1     
360            v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
361            v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
362            smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
363            oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
364            e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1 
365         ! Sample of adding more variables to apply lim_hdf---------
366         !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
367         !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
368         !-----------------------------------------------------------
369            DO jk = 1, nlay_i
370               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 
371            END DO
372         END DO
373
374         ato_i  (:,:) = zhdfptab(:,:,jm)
375
376         !------------------------------------------------------------------------------!
377         ! limit ice properties after transport                           
378         !------------------------------------------------------------------------------!
379!!gm & cr   :  MAX should not be active if adv scheme is positive !
380         DO jl = 1, jpl
381            DO jj = 1, jpj
382               DO ji = 1, jpi
383                  v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) )
384                  v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) )
385                  smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) )
386                  oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) )
387                  a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) )
388                  e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) )
389               END DO
390            END DO
391
392            DO jk = 1, nlay_i
393               DO jj = 1, jpj
394                  DO ji = 1, jpi
395                     e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) )
396                  END DO
397               END DO
398            END DO
399         END DO
400!!gm & cr
401
402         ! --- diags ---
403         DO jj = 1, jpj
404            DO ji = 1, jpi
405               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice
406               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice
407
408               diag_trp_vi (ji,jj) = SUM(   v_i(ji,jj,:) -  zviold(ji,jj,:) ) * r1_rdtice
409               diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice
410               diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice
411            END DO
412         END DO
413
414         ! zap small areas
415         CALL lim_var_zapsmall
416
417         !--- Thickness correction in case too high --------------------------------------------------------
418         DO jl = 1, jpl
419            DO jj = 1, jpj
420               DO ji = 1, jpi
421
422                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
423
424                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
425                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
426                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
427                     
428                     zvi  = v_i  (ji,jj,jl)
429                     zvs  = v_s  (ji,jj,jl)
430                     zsmv = smv_i(ji,jj,jl)
431                     zes  = e_s  (ji,jj,1,jl)
432                     zei  = SUM( e_i(ji,jj,1:nlay_i,jl) )
433
434                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 
435
436                     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. &
437                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
438
439                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
440                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
441
442                        ! small correction due to *rswitch for a_i
443                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
444                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
445                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
446                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
447                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
448
449                        ! Update mass fluxes
450                        wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
451                        wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
452                        sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
453                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
454                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
455
456                     ENDIF
457
458                  ENDIF
459
460               END DO
461            END DO
462         END DO
463         ! -------------------------------------------------
464         
465         !--------------------------------------
466         ! Impose a_i < amax in mono-category
467         !--------------------------------------
468         !
469         IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
470            DO jj = 1, jpj
471               DO ji = 1, jpi
472                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax_2d(ji,jj) )
473               END DO
474            END DO
475         ENDIF
476
477         ! --- agglomerate variables -----------------
478         vt_i (:,:) = 0._wp
479         vt_s (:,:) = 0._wp
480         at_i (:,:) = 0._wp
481         DO jl = 1, jpl
482            DO jj = 1, jpj
483               DO ji = 1, jpi
484                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl)
485                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl)
486                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
487               END DO
488            END DO
489         END DO
490
491         ! --- open water = 1 if at_i=0 --------------------------------
492         DO jj = 1, jpj
493            DO ji = 1, jpi
494               rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
495               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
496            END DO
497         END DO     
498
499         ! conservation test
500         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
501
502      ENDIF
503
504      ! -------------------------------------------------
505      ! control prints
506      ! -------------------------------------------------
507      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
508      !
509      CALL wrk_dealloc( jpi,jpj,            zsm, zatold, zeiold, zesold )
510      CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
511      CALL wrk_dealloc( jpi,jpj,1,          z0opw )
512      CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
513      CALL wrk_dealloc( jpi,jpj,jpl,        zviold, zvsold, zhimax, zsmvold )
514      CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars+nlay_i)+1,zhdfptab)
515      !
516      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
517
518   END SUBROUTINE lim_trp
519
520#else
521   !!----------------------------------------------------------------------
522   !!   Default option         Empty Module                No sea-ice model
523   !!----------------------------------------------------------------------
524CONTAINS
525   SUBROUTINE lim_trp        ! Empty routine
526   END SUBROUTINE lim_trp
527#endif
528   !!======================================================================
529END MODULE limtrp
530
Note: See TracBrowser for help on using the repository browser.