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

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 3963

Last change on this file since 3963 was 3963, checked in by clem, 11 years ago

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116
  • Property svn:keywords set to Id
File size: 33.4 KB
Line 
1MODULE limtrp
2   !!======================================================================
3   !!                       ***  MODULE limtrp   ***
4   !! LIM transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6   !! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code
7   !!            3.0  ! 2005-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_trp      : advection/diffusion process of sea ice
15   !!----------------------------------------------------------------------
16   USE phycst          ! physical constant
17   USE dom_oce         ! ocean domain
18   USE sbc_oce         ! ocean surface boundary condition
19   USE par_ice         ! LIM-3 parameter
20   USE dom_ice         ! LIM-3 domain
21   USE ice             ! LIM-3 variables
22   USE limadv          ! LIM-3 advection
23   USE limhdf          ! LIM-3 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     ! to use key_nosignedzero
30   USE limvar          ! clem for ice thickness correction
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   lim_trp    ! called by ice_step
36
37   REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values
38   REAL(wp)  ::   epsi03 = 1.e-03_wp 
39   REAL(wp)  ::   epsi10 = 1.e-10_wp 
40   REAL(wp)  ::   epsi16 = 1.e-16_wp
41   REAL(wp)  ::   epsi20 = 1.e-20_wp
42   REAL(wp)  ::   rzero  = 0._wp   
43   REAL(wp)  ::   rone   = 1._wp
44
45   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) ::   zs0e
46
47   !! * Substitution
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE lim_trp( kt ) 
57      !!-------------------------------------------------------------------
58      !!                   ***  ROUTINE lim_trp ***
59      !!                   
60      !! ** purpose : advection/diffusion process of sea ice
61      !!
62      !! ** method  : variables included in the process are scalar,   
63      !!     other values are considered as second order.
64      !!     For advection, a second order Prather scheme is used. 
65      !!
66      !! ** action :
67      !!---------------------------------------------------------------------
68      INTEGER, INTENT(in) ::   kt   ! number of iteration
69      !
70      INTEGER  ::   ji, jj, jk, jl, jm, layer   ! dummy loop indices
71      INTEGER  ::   jbnd1, jbnd2
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 , zrtt          !   -      -
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 ! 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
88      REAL(wp) :: zdv, zda, zvi, zvs, zsmv
89      !!---------------------------------------------------------------------
90
91      CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
92      CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
93      CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e )
94
95      CALL wrk_alloc( jpi,jpj,jpl,zviold )   ! clem
96      CALL wrk_alloc( jpi,jpj,jpl,zaiold, zhimax )   ! clem
97
98      ! -------------------------------
99      !- check conservation (C Rousset)
100      IF (ln_limdiahsb) THEN
101         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
102         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
103         zchk_fw_b  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) )
104         zchk_fs_b  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) )
105      ENDIF
106      !- check conservation (C Rousset)
107      ! -------------------------------
108
109      IF( numit == nstart .AND. lwp ) THEN
110         WRITE(numout,*)
111         IF( ln_limdyn ) THEN   ;   WRITE(numout,*) 'lim_trp : Ice transport '
112         ELSE                   ;   WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
113         ENDIF
114         WRITE(numout,*) '~~~~~~~~~~~~'
115      ENDIF
116     
117      zsm(:,:) = area(:,:)
118
119      !                             !-------------------------------------!
120      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
121         !                          !-------------------------------------!
122         ! mass and salt flux init (clem)
123         zviold(:,:,:)  = v_i(:,:,:)
124
125         !--- Thickness correction init. (clem) -------------------------------
126         CALL lim_var_glo2eqv
127         zaiold(:,:,:) = a_i(:,:,:)
128         !---------------------------------------------------------------------
129         ! Record max of the surrounding ice thicknesses for correction in limupdate
130         ! in case advection creates ice too thick.
131         !---------------------------------------------------------------------
132         zhimax(:,:,:) = ht_i(:,:,:)
133         DO jl = 1, jpl
134            DO jj = 2, jpjm1
135               DO ji = 2, jpim1
136                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
137               END DO
138            END DO
139            CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
140         END DO
141         
142         !-------------------------
143         ! transported fields                                       
144         !-------------------------
145         ! Snow vol, ice vol, salt and age contents, area
146         zs0ow(:,:) = ato_i(:,:) * area(:,:)               ! Open water area
147         DO jl = 1, jpl
148            zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume
149            zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume
150            zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area
151            zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content
152            zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content
153            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content
154            zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content
155         END DO
156
157         !--------------------------
158         ! Advection of Ice fields  (Prather scheme)                                           
159         !--------------------------
160         ! If ice drift field is too fast, use an appropriate time step for advection.         
161         ! CFL test for stability
162         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )
163         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
164         IF(lk_mpp )   CALL mpp_max( zcfl )
165!!gm more readability:
166!         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
167!         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
168!         ENDIF
169!!gm end
170         initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
171         zusnit = 1.0 / REAL( initad ) 
172         IF( zcfl > 0.5 .AND. lwp )   &
173            WRITE(numout,*) 'lim_trp   : CFL violation at day ', nday, ', cfl = ', zcfl,   &
174               &                        ': the ice time stepping is split in two'
175
176         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
177            DO jk = 1,initad
178               CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
179                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
180               CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   &
181                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
182               DO jl = 1, jpl
183                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
184                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
185                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
186                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
187                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
188                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
189                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
190                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
191                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
192                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
193                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
194                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
195                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---     
196                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
197                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
198                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
199                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
200                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
201                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
202                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
203                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
204                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
205                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
206                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
207                  DO layer = 1, nlay_i                                                           !--- ice heat contents ---
208                     CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
209                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
210                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
211                     CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
212                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
213                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
214                  END DO
215               END DO
216            END DO
217         ELSE
218            DO jk = 1, initad
219               CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
220                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
221               CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0ow (:,:), sxopw(:,:),   &
222                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
223               DO jl = 1, jpl
224                  CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
225                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
226                  CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
227                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
228                  CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
229                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
230                  CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
231                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
232                  CALL lim_adv_y( zusnit, v_ice, rone, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
233                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
234                  CALL lim_adv_x( zusnit, u_ice, rzero , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
235                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
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         !  Remultiply everything by ice area
327         !-----------------------------------------
328         !clem zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) )
329         !clem DO jl = 1, jpl
330         !clem    zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) )    !!bug:  est-ce utile
331         !clem    zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) )    !!bug:  cf /area  juste apres
332         !clem    zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) )    !!bug:  cf /area  juste apres
333         !clem    zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) )
334         !clem    zs0a  (:,:,jl) = MAX( rzero, zs0a  (:,:,jl) * area(:,:) )    !! suppress both change le resultat
335         !clem    zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) )
336         !clem    DO jk = 1, nlay_i
337         !clem       zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) )
338         !clem    END DO ! jk
339         !clem END DO ! jl
340
341         !------------------------------------------------------------------------------!
342         ! 5) Update and limit ice properties after transport                           
343         !------------------------------------------------------------------------------!
344
345         !--------------------------------------------------
346         ! 5.1) Recover mean values over the grid squares.
347         !--------------------------------------------------
348
349         !clem DO jl = 1, jpl
350         !clem    DO jk = 1, nlay_i
351         !clem       DO jj = 1, jpj
352         !clem          DO ji = 1, jpi
353         !clem             zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) )
354         !clem          END DO
355         !clem       END DO
356         !clem    END DO
357         !clem END DO
358
359         !clem DO jj = 1, jpj
360         !clem    DO ji = 1, jpi
361         !clem       zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) )
362         !clem    END DO
363         !clem END DO
364
365         zs0at(:,:) = 0._wp
366         DO jl = 1, jpl
367            DO jj = 1, jpj
368               DO ji = 1, jpi
369                 !clem  zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl)/area(ji,jj) )
370                 !clem  zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl)/area(ji,jj) )
371                 !clem  zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl)/area(ji,jj) )
372                 !clem  zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl)/area(ji,jj) )
373                 !clem  zs0a  (ji,jj,jl) = MAX( rzero, zs0a  (ji,jj,jl)/area(ji,jj) )
374                 !clem  zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl)/area(ji,jj) )
375                  zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) )
376                  zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) )
377                  zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) )
378                  zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) )
379                  zs0a  (ji,jj,jl) = MAX( rzero, zs0a  (ji,jj,jl) )
380                  zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) )
381                  zs0at (ji,jj)    = zs0at(ji,jj) + zs0a(ji,jj,jl)
382               END DO
383            END DO
384         END DO
385
386         !---------------------------------------------------------
387         ! 5.2) Snow thickness, Ice thickness, Ice concentrations
388         !---------------------------------------------------------
389         DO jj = 1, jpj
390            DO ji = 1, jpi
391               zindb        = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - epsi10) )
392               zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp )
393               ato_i(ji,jj) = zs0ow(ji,jj)
394            END DO
395         END DO
396
397         !
398         !
399         DO jl = 1, jpl         ! Remove very small areas
400            DO jj = 1, jpj
401               DO ji = 1, jpi
402                  zvi = zs0ice(ji,jj,jl)
403                  zvs = zs0sn(ji,jj,jl)
404                  !
405                  zindb         = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - epsi10) )
406                  !
407                  !zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 )
408                  v_s(ji,jj,jl)  = zindb * zs0sn (ji,jj,jl) 
409                  v_i(ji,jj,jl)  = zindb * zs0ice(ji,jj,jl)
410                  !
411                  zindsn         = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
412                  zindic         = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
413                  zindb          = MAX( zindsn, zindic )
414                  !
415                  zs0a(ji,jj,jl) = zindb  * zs0a(ji,jj,jl) !ice concentration
416                  a_i (ji,jj,jl) = zs0a(ji,jj,jl)
417                  v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl)
418                  v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl)
419                  !
420                  ! Update mass fluxes (clem)
421                  rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 
422                  rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 
423              END DO
424            END DO
425         END DO
426
427         !--- Thickness correction in case too high (clem) --------------------------------------------------------
428         CALL lim_var_glo2eqv
429         DO jl = 1, jpl
430            DO jj = 1, jpj
431               DO ji = 1, jpi
432
433                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
434                     zvi = v_i(ji,jj,jl)
435                     zvs = v_s(ji,jj,jl)
436                     zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl)   
437                     !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl)   
438                     
439                     IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. &
440                        & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                         
441                        ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
442                        zindh   =  MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )
443                        a_i(ji,jj,jl)  = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 )
444                     ELSE
445                        ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) )
446                        zindh   =  MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )
447                        a_i(ji,jj,jl)  = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 )
448                     ENDIF
449
450                     ! small correction due to *zindh for a_i
451                     v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl)
452                     v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl)
453
454                     ! Update mass fluxes
455                     rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic
456                     rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn
457
458                  ENDIF
459
460                  diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice
461
462               END DO
463            END DO
464         END DO
465
466         ! ---
467         DO jj = 1, jpj
468            DO ji = 1, jpi
469         !      ato_i(ji,jj) = 1._wp - SUM( a_i(ji,jj,1:jpl) ) !clem@rm-ow-advection
470               zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless??
471            END DO
472         END DO
473
474         !----------------------
475         ! 5.3) Ice properties
476         !----------------------
477
478         zbigval = 1.e+13
479
480         DO jl = 1, jpl
481            DO jj = 1, jpj
482               DO ji = 1, jpi
483                  zsmv = zs0sm(ji,jj,jl)
484
485                  ! Switches and dummy variables
486                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 )
487                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 )
488                  zrtt            = 173.15 * rone 
489                  zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
490                  zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
491                  zindb           = MAX( zindsn, zindic )
492
493                  ! Ice salinity and age
494                  zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , &
495                     zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl)
496                  IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) smv_i(ji,jj,jl) = zindic*zsal
497
498                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / & 
499                     MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * a_i(ji,jj,jl)
500                  oa_i (ji,jj,jl)  = zindic*zage 
501
502                  ! Snow heat content
503                  ze              =  MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval )
504                  e_s(ji,jj,1,jl) = zindsn * ze     
505
506                  ! Update salt fluxes (clem)
507                  fsalt_res(ji,jj) = fsalt_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic / rdt_ice 
508               END DO !ji
509            END DO !jj
510         END DO ! jl
511
512         DO jl = 1, jpl
513            DO jk = 1, nlay_i
514               DO jj = 1, jpj
515                  DO ji = 1, jpi
516                     ! Ice heat content
517                     zindic          =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
518                     ze              =  MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval )
519                     e_i(ji,jj,jk,jl) = zindic * ze
520                  END DO !ji
521               END DO ! jj
522            END DO ! jk
523         END DO ! jl
524
525
526      ! --- agglomerate variables (clem) -----------------
527      vt_i (:,:) = 0._wp
528      vt_s (:,:) = 0._wp
529      at_i (:,:) = 0._wp
530      !
531      DO jl = 1, jpl
532         DO jj = 1, jpj
533            DO ji = 1, jpi
534               !
535               vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
536               vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
537               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
538               !
539               zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi16 ) )
540               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness
541            END DO
542         END DO
543      END DO
544      ! -------------------------------------------------
545
546
547
548      ENDIF
549
550      IF(ln_ctl) THEN   ! Control print
551         CALL prt_ctl_info(' ')
552         CALL prt_ctl_info(' - Cell values : ')
553         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
554         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :')
555         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :')
556         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :')
557         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :')
558         DO jl = 1, jpl
559            CALL prt_ctl_info(' ')
560            CALL prt_ctl_info(' - Category : ', ivar1=jl)
561            CALL prt_ctl_info('   ~~~~~~~~~~')
562            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ')
563            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ')
564            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ')
565            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ')
566            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ')
567            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ')
568            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ')
569            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ')
570            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ')
571            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ')
572            DO jk = 1, nlay_i
573               CALL prt_ctl_info(' ')
574               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
575               CALL prt_ctl_info('   ~~~~~~~')
576               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ')
577               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ')
578            END DO
579         END DO
580      ENDIF
581      ! -------------------------------
582      !- check conservation (C Rousset)
583      IF (ln_limdiahsb) THEN
584         zchk_fs  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
585         zchk_fw  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
586 
587         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
588         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
589
590         zchk_vmin = glob_min(v_i)
591         zchk_amax = glob_max(SUM(a_i,dim=3))
592         zchk_amin = glob_min(a_i)
593
594         IF(lwp) THEN
595            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limtrp) = ',(zchk_v_i * 86400.)
596            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * 86400.)
597            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(zchk_vmin * 1.e-3)
598            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',zchk_amin
599         ENDIF
600      ENDIF
601      !- check conservation (C Rousset)
602      ! -------------------------------
603      !
604      CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
605      CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
606      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )
607
608      CALL wrk_dealloc( jpi,jpj,jpl,zaiold, zhimax )   ! clem
609      !
610   END SUBROUTINE lim_trp
611
612#else
613   !!----------------------------------------------------------------------
614   !!   Default option         Empty Module                No sea-ice model
615   !!----------------------------------------------------------------------
616CONTAINS
617   SUBROUTINE lim_trp        ! Empty routine
618   END SUBROUTINE lim_trp
619#endif
620
621   !!======================================================================
622END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.