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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.

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