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

source: branches/NERC/dev_r5107_NOC_MEDUSA/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 5715

Last change on this file since 5715 was 5715, checked in by acc, 9 years ago

Branch NERC/dev_r5107_NOC_MEDUSA. Complete reset of svn keyword properties in a desperate attempt to make fcm_make behave.

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