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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 10 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config

  • Property svn:keywords set to Id
File size: 31.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 par_ice        ! ice parameter
20   USE dom_ice        ! ice domain
21   USE ice            ! ice variables
22   USE limadv         ! ice advection
23   USE limhdf         ! ice horizontal diffusion
24   USE in_out_manager ! I/O manager
25   USE lbclnk         ! lateral boundary conditions -- MPP exchanges
26   USE lib_mpp        ! MPP library
27   USE wrk_nemo       ! work arrays
28   USE prtctl         ! Print control
29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30   USE limvar          ! clem for ice thickness correction
31   USE timing          ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   lim_trp    ! called by ice_step
37
38   REAL(wp)  ::   epsi10 = 1.e-10_wp 
39   REAL(wp)  ::   rzero  = 0._wp   
40   REAL(wp)  ::   rone   = 1._wp
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, layer   ! dummy loop indices
66      INTEGER  ::   initad                  ! number of sub-timestep for the advection
67      INTEGER  ::   ierr                    ! error status
68      REAL(wp) ::   zindb  , zindsn , zindic, zindh, zinda      ! local scalar
69      REAL(wp) ::   zusvosn, zusvoic, zbigval     !   -      -
70      REAL(wp) ::   zcfl , zusnit                 !   -      -
71      REAL(wp) ::   ze   , zsal   , zage          !   -      -
72      !
73      REAL(wp), POINTER, DIMENSION(:,:)      ::   zui_u, zvi_v, zsm, zs0at, zs0ow
74      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi
75      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e
76      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
77      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax, zchk_umax ! Check errors (C Rousset)
78      ! mass and salt flux (clem)
79      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold   ! old ice volume...
80      ! correct ice thickness (clem)
81      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zaiold, zhimax   ! old ice concentration and thickness
82      REAL(wp) :: zdv, zda, zvi, zvs, zsmv
83      !!---------------------------------------------------------------------
84      IF( nn_timing == 1 )  CALL timing_start('limtrp')
85
86      CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
87      CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
88      CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e )
89
90      CALL wrk_alloc( jpi,jpj,jpl,zviold )   ! clem
91      CALL wrk_alloc( jpi,jpj,jpl,zaiold, zhimax )   ! clem
92
93      ! -------------------------------
94      !- check conservation (C Rousset)
95      IF( ln_limdiahsb ) THEN
96         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
97         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
98         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
99         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
100      ENDIF
101      !- check conservation (C Rousset)
102      ! -------------------------------
103
104      IF( numit == nstart .AND. lwp ) THEN
105         WRITE(numout,*)
106         IF( ln_limdyn ) THEN   ;   WRITE(numout,*) 'lim_trp : Ice transport '
107         ELSE                   ;   WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
108         ENDIF
109         WRITE(numout,*) '~~~~~~~~~~~~'
110      ENDIF
111     
112      zsm(:,:) = area(:,:)
113
114      !                             !-------------------------------------!
115      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
116         !                          !-------------------------------------!
117         ! mass and salt flux init (clem)
118         zviold(:,:,:)  = v_i(:,:,:)
119
120         !--- Thickness correction init. (clem) -------------------------------
121         CALL lim_var_glo2eqv
122         zaiold(:,:,:) = a_i(:,:,:)
123         !---------------------------------------------------------------------
124         ! Record max of the surrounding ice thicknesses for correction in limupdate
125         ! in case advection creates ice too thick.
126         !---------------------------------------------------------------------
127         zhimax(:,:,:) = ht_i(:,:,:)
128         DO jl = 1, jpl
129            DO jj = 2, jpjm1
130               DO ji = 2, jpim1
131                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
132                  !zhimax(ji,jj,jl) = ( ht_i(ji  ,jj  ,jl) * tmask(ji,  jj  ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) &
133                  !     &             + ht_i(ji-1,jj  ,jl) * tmask(ji-1,jj  ,1) + ht_i(ji  ,jj-1,jl) * tmask(ji  ,jj-1,1) &
134                  !     &             + ht_i(ji+1,jj  ,jl) * tmask(ji+1,jj  ,1) + ht_i(ji  ,jj+1,jl) * tmask(ji  ,jj+1,1) &
135                  !     &             + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) )
136               END DO
137            END DO
138            CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
139         END DO
140         
141         !-------------------------
142         ! transported fields                                       
143         !-------------------------
144         ! Snow vol, ice vol, salt and age contents, area
145         zs0ow(:,:) = ato_i(:,:) * area(:,:)               ! Open water area
146         DO jl = 1, jpl
147            zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume
148            zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume
149            zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area
150            zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content
151            zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content
152            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content
153            zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content
154         END DO
155
156         !--------------------------
157         ! Advection of Ice fields  (Prather scheme)                                           
158         !--------------------------
159         ! If ice drift field is too fast, use an appropriate time step for advection.         
160         ! CFL test for stability
161         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )
162         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
163         IF(lk_mpp )   CALL mpp_max( zcfl )
164!!gm more readability:
165!         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
166!         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
167!         ENDIF
168!!gm end
169         initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
170         zusnit = 1.0 / REAL( initad ) 
171         IF( zcfl > 0.5 .AND. lwp )   &
172            WRITE(numout,*) 'lim_trp   : CFL violation at day ', nday, ', cfl = ', zcfl,   &
173               &                        ': the ice time stepping is split in two'
174
175         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
176            DO jk = 1,initad
177               CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
178                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
179               CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   &
180                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
181               DO jl = 1, jpl
182                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
183                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
184                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
185                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
186                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
187                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
188                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
189                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
190                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
191                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
192                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
193                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
194                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---     
195                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
196                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
197                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
198                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
199                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
200                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
201                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
202                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
203                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
204                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
205                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
206                  DO layer = 1, nlay_i                                                           !--- ice heat contents ---
207                     CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
208                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
209                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
210                     CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
211                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
212                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
213                  END DO
214               END DO
215            END DO
216         ELSE
217            DO jk = 1, initad
218               CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
219                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
220               CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   &
221                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
222               DO jl = 1, jpl
223                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
224                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
225                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
226                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
227                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
228                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
229                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
230                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
231                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
232                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
233                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
234                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
235
236                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
237                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
238                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
239                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
240                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
241                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
242                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &
243                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
244                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
245                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
246                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
247                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
248                  DO layer = 1, nlay_i                                                           !--- ice heat contents ---
249                     CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
250                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
251                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
252                     CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
253                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
254                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
255                  END DO
256               END DO
257            END DO
258         ENDIF
259
260         !-------------------------------------------
261         ! Recover the properties from their contents
262         !-------------------------------------------
263         zs0ow(:,:) = zs0ow(:,:) / area(:,:)
264         DO jl = 1, jpl
265            zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)
266            zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)
267            zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)
268            zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)
269            zs0a  (:,:,jl) = zs0a  (:,:,jl) / area(:,:)
270            zs0c0 (:,:,jl) = zs0c0 (:,:,jl) / area(:,:)
271            DO jk = 1, nlay_i
272               zs0e(:,:,jk,jl) = zs0e(:,:,jk,jl) / area(:,:)
273            END DO
274         END DO
275
276         !------------------------------------------------------------------------------!
277         ! 4) Diffusion of Ice fields                 
278         !------------------------------------------------------------------------------!
279
280         !--------------------------------
281         !  diffusion of open water area
282         !--------------------------------
283         zs0at(:,:) = zs0a(:,:,1)      ! total ice fraction
284         DO jl = 2, jpl
285            zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)
286         END DO
287         !
288         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
289         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
290            DO ji = 1 , fs_jpim1   ! vector opt.
291               pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji  ,jj) ) ) )   &
292                  &        * ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)
293               pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji,jj  ) ) ) )   &
294                  &        * ( 1._wp - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)
295            END DO
296         END DO
297         !
298         CALL lim_hdf( zs0ow (:,:) )   ! Diffusion
299
300         !------------------------------------
301         !  Diffusion of other ice variables
302         !------------------------------------
303         DO jl = 1, jpl
304         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
305            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
306               DO ji = 1 , fs_jpim1   ! vector opt.
307                  pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj,jl) ) ) )   &
308                     &        * ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
309                  pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji,jj  ,jl) ) ) )   &
310                     &        * ( 1._wp - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
311               END DO
312            END DO
313
314            CALL lim_hdf( zs0ice (:,:,jl) )
315            CALL lim_hdf( zs0sn  (:,:,jl) )
316            CALL lim_hdf( zs0sm  (:,:,jl) )
317            CALL lim_hdf( zs0oi  (:,:,jl) )
318            CALL lim_hdf( zs0a   (:,:,jl) )
319            CALL lim_hdf( zs0c0  (:,:,jl) )
320            DO jk = 1, nlay_i
321               CALL lim_hdf( zs0e (:,:,jk,jl) )
322            END DO
323         END DO
324
325         !------------------------------------------------------------------------------!
326         ! 5) Update and limit ice properties after transport                           
327         !------------------------------------------------------------------------------!
328
329         !--------------------------------------------------
330         ! 5.1) Recover mean values over the grid squares.
331         !--------------------------------------------------
332         zs0at(:,:) = 0._wp
333         DO jl = 1, jpl
334            DO jj = 1, jpj
335               DO ji = 1, jpi
336                  zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) )
337                  zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) )
338                  zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) )
339                  zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) )
340                  zs0a  (ji,jj,jl) = MAX( rzero, zs0a  (ji,jj,jl) )
341                  zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) )
342                  zs0at (ji,jj)    = zs0at(ji,jj) + zs0a(ji,jj,jl)
343               END DO
344            END DO
345         END DO
346
347         !---------------------------------------------------------
348         ! 5.2) Snow thickness, Ice thickness, Ice concentrations
349         !---------------------------------------------------------
350         DO jj = 1, jpj
351            DO ji = 1, jpi
352               zindb        = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - epsi10) )
353               zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp )
354               ato_i(ji,jj) = zs0ow(ji,jj)
355            END DO
356         END DO
357
358         DO jl = 1, jpl         ! Remove very small areas
359            DO jj = 1, jpj
360               DO ji = 1, jpi
361                  zvi = zs0ice(ji,jj,jl)
362                  zvs = zs0sn(ji,jj,jl)
363                  !
364                  zindb         = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - epsi10) )
365                  !
366                  v_s(ji,jj,jl)  = zindb * zs0sn (ji,jj,jl) 
367                  v_i(ji,jj,jl)  = zindb * zs0ice(ji,jj,jl)
368                  !
369                  zindsn         = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
370                  zindic         = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
371                  zindb          = MAX( zindsn, zindic )
372                  !
373                  zs0a(ji,jj,jl) = zindb  * zs0a(ji,jj,jl) !ice concentration
374                  a_i (ji,jj,jl) = zs0a(ji,jj,jl)
375                  v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl)
376                  v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl)
377                  !
378                  ! Update mass fluxes (clem)
379                  rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 
380                  rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 
381              END DO
382            END DO
383         END DO
384
385         !--- Thickness correction in case too high (clem) --------------------------------------------------------
386         CALL lim_var_glo2eqv
387         DO jl = 1, jpl
388            DO jj = 1, jpj
389               DO ji = 1, jpi
390
391                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
392                     zvi = v_i(ji,jj,jl)
393                     zvs = v_s(ji,jj,jl)
394                     zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl)   
395                     !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl)   
396                     
397                     zindh = 1._wp
398                     IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. &
399                        & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                         
400                        ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
401                        zindh   =  MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )
402                        a_i(ji,jj,jl)  = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 )
403                     ELSE
404                        ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) )
405                        zindh   =  MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )
406                        a_i(ji,jj,jl)  = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 )
407                     ENDIF
408
409                     ! small correction due to *zindh for a_i
410                     v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl)
411                     v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl)
412
413                     ! Update mass fluxes
414                     rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic
415                     rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn
416
417                  ENDIF
418
419                  diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice
420
421               END DO
422            END DO
423         END DO
424
425         ! ---
426         DO jj = 1, jpj
427            DO ji = 1, jpi
428               zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless??
429            END DO
430         END DO
431
432         !----------------------
433         ! 5.3) Ice properties
434         !----------------------
435
436         zbigval = 1.e+13
437
438         DO jl = 1, jpl
439            DO jj = 1, jpj
440               DO ji = 1, jpi
441                  zsmv = zs0sm(ji,jj,jl)
442
443                  ! Switches and dummy variables
444                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi10 )
445                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi10 )
446                  zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
447                  zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
448                  zindb           = MAX( zindsn, zindic )
449
450                  ! Ice salinity and age
451                  !clem zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj), zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl)
452                  IF(  num_sal == 2  ) THEN
453                     smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) )
454                  ENDIF
455
456                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) ), 0._wp  ) * a_i(ji,jj,jl)
457                  oa_i (ji,jj,jl)  = zindic * zage 
458
459                  ! Snow heat content
460                  ze              =  MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval )
461                  e_s(ji,jj,1,jl) = zindsn * ze     
462
463                  ! Update salt fluxes (clem)
464                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
465               END DO !ji
466            END DO !jj
467         END DO ! jl
468
469         DO jl = 1, jpl
470            DO jk = 1, nlay_i
471               DO jj = 1, jpj
472                  DO ji = 1, jpi
473                     ! Ice heat content
474                     zindic          =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
475                     ze              =  MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval )
476                     e_i(ji,jj,jk,jl) = zindic * ze
477                  END DO !ji
478               END DO ! jj
479            END DO ! jk
480         END DO ! jl
481
482
483      ! --- agglomerate variables (clem) -----------------
484      vt_i (:,:) = 0._wp
485      vt_s (:,:) = 0._wp
486      at_i (:,:) = 0._wp
487      !
488      DO jl = 1, jpl
489         DO jj = 1, jpj
490            DO ji = 1, jpi
491               !
492               vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
493               vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
494               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
495               !
496               zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi10 ) )
497               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness
498            END DO
499         END DO
500      END DO
501      ! -------------------------------------------------
502
503
504
505      ENDIF
506
507      IF(ln_ctl) THEN   ! Control print
508         CALL prt_ctl_info(' ')
509         CALL prt_ctl_info(' - Cell values : ')
510         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
511         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :')
512         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :')
513         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :')
514         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :')
515         DO jl = 1, jpl
516            CALL prt_ctl_info(' ')
517            CALL prt_ctl_info(' - Category : ', ivar1=jl)
518            CALL prt_ctl_info('   ~~~~~~~~~~')
519            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ')
520            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ')
521            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ')
522            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ')
523            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ')
524            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ')
525            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ')
526            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ')
527            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ')
528            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ')
529            DO jk = 1, nlay_i
530               CALL prt_ctl_info(' ')
531               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
532               CALL prt_ctl_info('   ~~~~~~~')
533               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ')
534               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ')
535            END DO
536         END DO
537      ENDIF
538      ! -------------------------------
539      !- check conservation (C Rousset)
540      IF( ln_limdiahsb ) THEN
541         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
542         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
543 
544         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
545         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
546
547         zchk_vmin = glob_min(v_i)
548         zchk_amax = glob_max(SUM(a_i,dim=3))
549         zchk_amin = glob_min(a_i)
550         zchk_umax = glob_max(SQRT(u_ice**2 + v_ice**2))
551
552         IF(lwp) THEN
553            IF ( ABS( zchk_v_i   ) >  1.e-5 ) THEN
554               WRITE(numout,*) 'violation volume [m3/day]     (limtrp) = ',(zchk_v_i * rday)
555               WRITE(numout,*) 'u_ice max [m/s]               (limtrp) = ',zchk_umax
556               WRITE(numout,*) 'number of time steps          (limtrp) =',kt
557            ENDIF
558            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday)
559            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(zchk_vmin * 1.e-3)
560            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',zchk_amin
561         ENDIF
562      ENDIF
563      !- check conservation (C Rousset)
564      ! -------------------------------
565      !
566      CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
567      CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
568      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )
569
570      CALL wrk_dealloc( jpi,jpj,jpl,zaiold, zhimax )   ! clem
571      !
572      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
573   END SUBROUTINE lim_trp
574
575#else
576   !!----------------------------------------------------------------------
577   !!   Default option         Empty Module                No sea-ice model
578   !!----------------------------------------------------------------------
579CONTAINS
580   SUBROUTINE lim_trp        ! Empty routine
581   END SUBROUTINE lim_trp
582#endif
583
584   !!======================================================================
585END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.