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

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 7319

Last change on this file since 7319 was 7319, checked in by vancop, 7 years ago

Sea ice model infrastructure for melt ponds part 2 (advection, transport in thickness space, ridging/rafting)

  • Property svn:keywords set to Id
File size: 34.1 KB
Line 
1MODULE limtrp
2   !!======================================================================
3   !!                       ***  MODULE limtrp   ***
4   !! LIM transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6   !! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code
7   !!            3.0  ! 2005-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_trp      : advection/diffusion process of sea ice
15   !!----------------------------------------------------------------------
16   USE phycst         ! physical constant
17   USE dom_oce        ! ocean domain
18   USE sbc_oce        ! ocean surface boundary condition
19   USE ice            ! ice variables
20   USE limhdf         ! ice horizontal diffusion
21   USE limvar         !
22   USE limadv_prather ! advection scheme (Prather)
23   USE limadv_umx     ! advection scheme (ultimate-macho)
24   !
25   USE in_out_manager ! I/O manager
26   USE lbclnk         ! lateral boundary conditions -- MPP exchanges
27   USE lib_mpp        ! MPP library
28   USE wrk_nemo       ! work arrays
29   USE prtctl         ! Print control
30   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
31   USE timing         ! Timing
32   USE limcons        ! conservation tests
33   USE limctl         ! control prints
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   lim_trp    ! called by sbcice_lim
39
40   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2 
41
42   !! * Substitution
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE lim_trp( kt ) 
52      !!-------------------------------------------------------------------
53      !!                   ***  ROUTINE lim_trp ***
54      !!                   
55      !! ** purpose : advection/diffusion process of sea ice
56      !!
57      !! ** method  : variables included in the process are scalar,   
58      !!     other values are considered as second order.
59      !!     For advection, one can choose between
60      !!     a) an Ultimate-Macho scheme (whose order is defined by nn_limadv_ord) => nn_limadv=0
61      !!     b) and a second order Prather scheme => nn_limadv=-1
62      !!
63      !! ** action :
64      !!---------------------------------------------------------------------
65      INTEGER, INTENT(in) ::   kt   ! number of iteration
66      !
67      INTEGER  ::   ji, jj, jk, jm, jl, jt  ! dummy loop indices
68      INTEGER  ::   initad                  ! number of sub-timestep for the advection
69      REAL(wp) ::   zcfl , zusnit           !   -      -
70      CHARACTER(len=80) :: cltmp
71      !
72      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
73      REAL(wp) ::    zdv, zda
74      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold, zsmvold 
75      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax, zviold, zvsold
76      ! --- diffusion --- !
77      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhdfptab
78      ! MV MP 2016
79      ! With melt ponds, we have to diffuse them
80      ! We hard code the number of variables to diffuse
81      ! since we can't put an IF ( ln_limMP ) for a declaration
82      ! ideally, the ihdf_vars should probably be passed as an argument and
83      ! defined somewhere depending on ln_limMP
84      ! END MV MP 2016
85      INTEGER , PARAMETER                    ::   ihdf_vars  = 8 ! Number of variables in which we apply horizontal diffusion
86                                                                 !  inside limtrp for each ice category , not counting the
87                                                                 !  variables corresponding to ice_layers
88      ! --- ultimate macho only --- !
89      REAL(wp)                               ::   zdt
90      REAL(wp), POINTER, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box
91      ! --- prather only --- !
92      REAL(wp), POINTER, DIMENSION(:,:)      ::   zarea
93      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
94      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
95      ! MV MP 2016
96      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ap , z0vp
97      REAL(wp) ::   za_old
98      ! END MV MP 2016
99      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
100      !!
101      !!---------------------------------------------------------------------
102      IF( nn_timing == 1 )  CALL timing_start('limtrp')
103
104      CALL wrk_alloc( jpi,jpj,                            zatold, zeiold, zesold, zsmvold )
105      CALL wrk_alloc( jpi,jpj,jpl,                        zhimax, zviold, zvsold )
106      CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1, zhdfptab)
107 
108      IF( kt == nit000 .AND. lwp ) THEN
109         WRITE(numout,*)''
110         WRITE(numout,*)'limtrp'
111         WRITE(numout,*)'~~~~~~'
112         ncfl = 0                ! nb of time step with CFL > 1/2
113      ENDIF
114     
115      CALL lim_var_agg( 1 ) ! integrated values + ato_i
116
117      !-------------------------------------!
118      !   Advection of sea ice properties   !
119      !-------------------------------------!
120
121      ! conservation test
122      IF( ln_limdiachk )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
123     
124      ! store old values for diag
125      zviold = v_i
126      zvsold = v_s
127      zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 )
128      zeiold (:,:) = et_i
129      zesold (:,:) = et_s 
130
131      !--- Thickness correction init. --- !
132      zatold(:,:) = at_i
133      DO jl = 1, jpl
134         DO jj = 1, jpj
135            DO ji = 1, jpi
136               rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
137               ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
138               ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
139            END DO
140         END DO
141      END DO
142      ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- !
143      zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
144      DO jl = 1, jpl
145         DO jj = 2, jpjm1
146            DO ji = 2, jpim1
147               zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )
148            END DO
149         END DO
150         CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
151      END DO
152         
153      ! --- If ice drift field is too fast, use an appropriate time step for advection --- !       
154      zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )    ! CFL test for stability
155      zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
156      IF( lk_mpp )   CALL mpp_max( zcfl )
157     
158      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
159      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
160      ENDIF
161     
162!!      IF( zcfl > 0.5_wp .AND. lwp ) THEN
163!!         ncfl = ncfl + 1
164!!         IF( ncfl > 0 ) THEN   
165!!            WRITE(cltmp,'(i6.1)') ncfl
166!!            CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
167!!         ENDIF
168!!      ENDIF
169
170      SELECT CASE ( nn_limadv )
171         
172                       !=============================!
173      CASE ( 0 )       !==  Ultimate-MACHO scheme  ==!                   
174                       !=============================!
175     
176         CALL wrk_alloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
177     
178         IF( kt == nit000 .AND. lwp ) THEN
179            WRITE(numout,*)''
180            WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme'
181            WRITE(numout,*)'~~~~~~~~~~~'
182         ENDIF
183         !
184         zdt = rdt_ice / REAL(initad)
185         
186         ! transport
187         zudy(:,:) = u_ice(:,:) * e2u(:,:)
188         zvdx(:,:) = v_ice(:,:) * e1v(:,:)
189         
190         ! define velocity for advection: u*grad(H)
191         DO jj = 2, jpjm1
192            DO ji = fs_2, fs_jpim1
193               IF    ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
194               ELSEIF( u_ice(ji,jj)                  >  0._wp ) THEN   ;   zcu_box(ji,jj) = u_ice(ji-1,jj)
195               ELSE                                                    ;   zcu_box(ji,jj) = u_ice(ji  ,jj)
196               ENDIF
197               
198               IF    ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
199               ELSEIF( v_ice(ji,jj)                  >  0._wp ) THEN   ;   zcv_box(ji,jj) = v_ice(ji,jj-1)
200               ELSE                                                    ;   zcv_box(ji,jj) = v_ice(ji,jj  )
201               ENDIF
202            END DO
203         END DO
204         
205         ! advection
206         DO jt = 1, initad
207            CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) )       ! Open water area
208            DO jl = 1, jpl
209               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) )      ! Ice area
210               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) )      ! Ice  volume
211               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) )    ! Salt content
212               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) )    ! Age content
213               DO jk = 1, nlay_i
214                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) )   ! Ice  heat content
215               END DO
216               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume
217               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content
218               ! MV MP 2016
219               IF ( ln_limMP ) THEN
220                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_ip(:,:,jl) )  ! Melt pond fraction
221                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_ip(:,:,jl) )  ! Melt pond volume
222               ENDIF
223               ! END MV MP 2016
224            END DO
225         END DO
226         !
227         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
228         DO jl = 2, jpl
229            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
230         END DO
231         !
232         CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
233         
234                       !=============================!
235      CASE ( -1 )      !==     Prather scheme      ==!                   
236                       !=============================!
237
238         CALL wrk_alloc( jpi,jpj,            zarea )
239         CALL wrk_alloc( jpi,jpj,1,          z0opw )
240         CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
241         CALL wrk_alloc( jpi,jpj,jpl,        z0ap , z0vp )
242         CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
243         
244         IF( kt == nit000 .AND. lwp ) THEN
245            WRITE(numout,*)''
246            WRITE(numout,*)'lim_adv_xy : Prather advection scheme'
247            WRITE(numout,*)'~~~~~~~~~~~'
248         ENDIF
249         
250         zarea(:,:) = e12t(:,:)
251         
252         !-------------------------
253         ! transported fields                                       
254         !-------------------------
255         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area
256         DO jl = 1, jpl
257            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
258            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
259            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
260            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
261            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
262            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
263            DO jk = 1, nlay_i
264               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
265            END DO
266            ! MV MP 2016
267            IF ( ln_limMP ) THEN
268               z0ap  (:,:,jl)  = a_ip (:,:,jl) * e12t(:,:) ! Melt pond fraction
269               z0vp  (:,:,jl)  = v_ip (:,:,jl) * e12t(:,:) ! Melt pond volume
270            ENDIF
271            ! END MV MP 2016
272
273         END DO
274         
275         
276         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
277            DO jt = 1, initad
278               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
279                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
280               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
281                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
282               DO jl = 1, jpl
283                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
284                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
285                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
286                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
287                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
288                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
289                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
290                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
291                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
292                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
293                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
294                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
295                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
296                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
297                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
298                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
299                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
300                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
301                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
302                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
303                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
304                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
305                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
306                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
307                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
308                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
309                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
310                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
311                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
312                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
313                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
314                  END DO
315                  ! MV MP 2016
316                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
317                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
318                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
319                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
320                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   --
321                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
322                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
323                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
324                  ! END MV MP 2016
325               END DO
326            END DO
327         ELSE
328            DO jt = 1, initad
329               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
330                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
331               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
332                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
333               DO jl = 1, jpl
334                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
335                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
336                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
337                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
338                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
339                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
340                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
341                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
342                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
343                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
344                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
345                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
346                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
347                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
348                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
349                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
350                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
351                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
352                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
353                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
354                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
355                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
356                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
357                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
358                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
359                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
360                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
361                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
362                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
363                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
364                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
365                  END DO
366                  ! MV MP 2016
367                  IF ( ln_limMP ) THEN
368                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
369                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
370                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &
371                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
372                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   ---
373                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
374                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &
375                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
376                  ENDIF
377                  ! END MV MP 2016
378               END DO
379            END DO
380         ENDIF
381         
382         !-------------------------------------------
383         ! Recover the properties from their contents
384         !-------------------------------------------
385         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
386         DO jl = 1, jpl
387            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
388            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
389            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
390            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
391            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
392            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
393            DO jk = 1, nlay_i
394               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
395            END DO
396            ! MV MP 2016
397            IF ( ln_limMP ) THEN
398               a_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e12t(:,:)
399               v_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e12t(:,:)
400            ENDIF
401            ! END MV MP 2016
402         END DO
403         
404         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
405         DO jl = 2, jpl
406            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
407         END DO
408         
409         CALL wrk_dealloc( jpi,jpj,            zarea )
410         CALL wrk_dealloc( jpi,jpj,1,          z0opw )
411         CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
412         ! MV MP 2016
413         CALL wrk_dealloc( jpi,jpj,jpl,        z0ap , z0vp )
414         ! END MV MP 2016
415         CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
416
417      END SELECT
418     
419      !------------------------------!
420      ! Diffusion of Ice fields                 
421      !------------------------------!
422      IF( nn_ahi0 /= -1 .AND. nn_limdyn == 2 ) THEN
423         !
424         ! --- Prepare diffusion for variables with categories --- !
425         !     mask eddy diffusivity coefficient at ocean U- and V-points
426         jm=1
427         DO jl = 1, jpl
428            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
429               DO ji = 1 , fs_jpim1
430                  pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,  jl ) ) ) )   &
431                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj)
432                  pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,  jj,  jl ) ) ) )   &
433                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,  jj+1,jl ) ) ) ) * ahiv(ji,jj)
434               END DO
435            END DO
436
437            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1
438            zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1
439            zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1
440            zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1
441            zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
442            zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
443            ! MV MP 2016
444            IF ( ln_limMP ) THEN
445               zhdfptab(:,:,jm)= a_ip  (:,:,  jl); jm = jm + 1
446               zhdfptab(:,:,jm)= v_ip  (:,:,  jl); jm = jm + 1
447            ENDIF
448            ! END MV MP 2016
449            ! Sample of adding more variables to apply lim_hdf (ihdf_vars must be increased)
450            !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1 
451            !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1
452            DO jk = 1, nlay_i
453              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
454            END DO
455         END DO
456
457         ! --- Prepare diffusion for open water area --- !
458         !     mask eddy diffusivity coefficient at ocean U- and V-points
459         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
460            DO ji = 1 , fs_jpim1
461               pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
462                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
463               pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
464                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
465            END DO
466         END DO
467         !
468         zhdfptab(:,:,jm)= ato_i  (:,:);
469
470         ! --- Apply diffusion --- !
471         CALL lim_hdf( zhdfptab, ihdf_vars )
472
473         ! --- Recover properties --- !
474         jm=1
475         DO jl = 1, jpl
476            a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
477            v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
478            v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
479            smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
480            oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
481            e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1
482            ! MV MP 2016
483            IF ( ln_limMP ) THEN
484               a_ip (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
485               v_ip (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
486            ENDIF
487            ! Sample of adding more variables to apply lim_hdf
488            !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
489            !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
490            DO jk = 1, nlay_i
491               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1
492            END DO
493         END DO
494         ato_i  (:,:) = zhdfptab(:,:,jm)
495             
496      ENDIF
497
498      ! --- diags ---
499      DO jj = 1, jpj
500         DO ji = 1, jpi
501            diag_trp_ei (ji,jj) = ( SUM( e_i  (ji,jj,1:nlay_i,:) ) -  zeiold(ji,jj) ) * r1_rdtice
502            diag_trp_es (ji,jj) = ( SUM( e_s  (ji,jj,1:nlay_s,:) ) -  zesold(ji,jj) ) * r1_rdtice
503            diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:)          ) - zsmvold(ji,jj) ) * r1_rdtice
504            diag_trp_vi (ji,jj) =   SUM(   v_i(ji,jj,:)            -  zviold(ji,jj,:) ) * r1_rdtice
505            diag_trp_vs (ji,jj) =   SUM(   v_s(ji,jj,:)            -  zvsold(ji,jj,:) ) * r1_rdtice
506         END DO
507      END DO
508     
509      IF( nn_limdyn == 2) THEN
510
511         ! zap small areas
512         CALL lim_var_zapsmall
513           
514         !--- Thickness correction in case too high --- !
515         DO jl = 1, jpl
516            DO jj = 1, jpj
517               DO ji = 1, jpi
518                 
519                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
520                     
521                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
522                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
523                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
524                     
525                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 
526                     
527                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
528                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
529                       
530                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
531                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
532                       
533                        ! small correction due to *rswitch for a_i
534                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
535                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
536                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
537                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
538                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
539
540                        ! MV MP 2016
541                        IF ( ln_limMP ) THEN
542                           a_ip (ji,jj,jl)        = rswitch * a_ip (ji,jj,jl)
543                           v_ip (ji,jj,jl)        = rswitch * v_ip (ji,jj,jl)
544                        ENDIF
545                        ! END MV MP 2016
546                                               
547                     ENDIF
548                     
549                  ENDIF
550               
551               END DO
552            END DO
553         END DO
554         
555         ! Force the upper limit of ht_i to always be < hi_max (99 m).
556         DO jj = 1, jpj
557            DO ji = 1, jpi
558               ! MV MP 2016
559               za_old = a_i(ji,jj,jpl)
560               ! END MV MP 2016
561               rswitch         = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
562               ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
563               a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
564               ! MV MP 2016
565               IF ( ln_limMP ) THEN
566                  ! correct pond fraction to avoid a_ip > a_i
567                  a_ip(ji,jj,jpl) = a_ip(ji,jj,jpl) * a_i(ji,jj,jpl) / MAX( za_old , epsi20 ) * rswitch
568               ENDIF
569               ! END MP 2016
570            END DO
571         END DO
572
573
574      ENDIF
575         
576      !------------------------------------------------------------
577      ! Impose a_i < amax if no ridging/rafting or in mono-category
578      !------------------------------------------------------------
579      !
580      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
581      IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2
582         DO jl = 1, jpl
583            DO jj = 1, jpj
584               DO ji = 1, jpi
585                  rswitch       = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) )
586                  zda           = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp )  &
587                     &                    * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 )
588                  a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda
589               END DO
590            END DO
591         END DO
592      ENDIF
593     
594      ! --- agglomerate variables -----------------
595      vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
596      vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
597      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
598
599      ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i)
600      IF ( ln_limMP ) THEN
601          at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 )
602          vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 )
603      ENDIF
604      ! END MP 2016
605     
606      ! --- open water = 1 if at_i=0 --------------------------------
607      WHERE( at_i == 0._wp ) ato_i = 1._wp 
608     
609      ! conservation test
610      IF( ln_limdiachk )   CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
611       
612      ! -------------------------------------------------
613      ! control prints
614      ! -------------------------------------------------
615      IF( ln_limctl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
616      !
617      CALL wrk_dealloc( jpi,jpj,                            zatold, zeiold, zesold, zsmvold )
618      CALL wrk_dealloc( jpi,jpj,jpl,                        zhimax, zviold, zvsold )
619      CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1, zhdfptab)
620      !
621      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
622      !
623   END SUBROUTINE lim_trp
624
625#else
626   !!----------------------------------------------------------------------
627   !!   Default option         Empty Module                No sea-ice model
628   !!----------------------------------------------------------------------
629CONTAINS
630   SUBROUTINE lim_trp        ! Empty routine
631   END SUBROUTINE lim_trp
632#endif
633
634   !!======================================================================
635END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.