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

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 4220

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

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