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 @ 4635

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

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