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/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 5059

Last change on this file since 5059 was 5059, checked in by clem, 9 years ago

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