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.
iceadv.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90 @ 8411

Last change on this file since 8411 was 8411, checked in by clem, 7 years ago

continue changing names

File size: 29.2 KB
Line 
1MODULE iceadv
2   !!======================================================================
3   !!                       ***  MODULE iceadv   ***
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   !!   ice_adv      : advection/diffusion process of sea ice
15   !!----------------------------------------------------------------------
16   USE phycst         ! physical constant
17   USE dom_oce        ! ocean domain
18   USE sbc_oce , ONLY : nn_fsbc
19   USE ice            ! ice variables
20   USE limvar         !
21   USE iceadv_prather ! advection scheme (Prather)
22   USE iceadv_umx     ! advection scheme (ultimate-macho)
23   !
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 timing         ! Timing
31   USE icecons        ! conservation tests
32   USE icectl         ! control prints
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   ice_adv    ! called by icestp
38
39   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2 
40
41   !! * Substitution
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
45   !! $Id: iceadv.F90 8373 2017-07-25 17:44:54Z clem $
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE ice_adv( kt ) 
51      !!-------------------------------------------------------------------
52      !!                   ***  ROUTINE ice_adv ***
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, one can choose between
59      !!     a) an Ultimate-Macho scheme (whose order is defined by nn_limadv_ord) => nn_limadv=0
60      !!     b) and a second order Prather scheme => nn_limadv=-1
61      !!
62      !! ** action :
63      !!---------------------------------------------------------------------
64      INTEGER, INTENT(in) ::   kt   ! number of iteration
65      !
66      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
67      INTEGER  ::   initad                  ! number of sub-timestep for the advection
68      REAL(wp) ::   zcfl , zusnit           !   -      -
69      CHARACTER(len=80) :: cltmp
70      !
71      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
72      REAL(wp) ::    zdv, zda
73      REAL(wp), DIMENSION(jpi,jpj)           ::   zatold, zeiold, zesold, zsmvold 
74      REAL(wp), DIMENSION(jpi,jpj,jpl)       ::   zhimax, zviold, zvsold
75      ! --- ultimate macho only --- !
76      REAL(wp)                               ::   zdt
77      REAL(wp), POINTER, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box
78      ! --- prather only --- !
79      REAL(wp), POINTER, DIMENSION(:,:)      ::   zarea
80      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
81      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
82      ! MV MP 2016
83      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ap , z0vp
84      REAL(wp) ::   za_old
85      ! END MV MP 2016
86      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
87      !!
88      !!---------------------------------------------------------------------
89      IF( nn_timing == 1 )  CALL timing_start('iceadv')
90 
91      IF( kt == nit000 .AND. lwp ) THEN
92         WRITE(numout,*)''
93         WRITE(numout,*)'iceadv'
94         WRITE(numout,*)'~~~~~~'
95         ncfl = 0                ! nb of time step with CFL > 1/2
96      ENDIF
97     
98      CALL lim_var_agg( 1 ) ! integrated values + ato_i
99
100      !-------------------------------------!
101      !   Advection of sea ice properties   !
102      !-------------------------------------!
103
104      ! conservation test
105      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
106     
107      ! store old values for diag
108      zviold = v_i
109      zvsold = v_s
110      zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 )
111      zeiold (:,:) = et_i
112      zesold (:,:) = et_s 
113
114      !--- Thickness correction init. --- !
115      zatold(:,:) = at_i
116      DO jl = 1, jpl
117         DO jj = 1, jpj
118            DO ji = 1, jpi
119               rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
120               ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
121               ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
122            END DO
123         END DO
124      END DO
125      ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- !
126      zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
127      DO jl = 1, jpl
128         DO jj = 2, jpjm1
129            DO ji = 2, jpim1
130               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) )
131            END DO
132         END DO
133      END DO
134      CALL lbc_lnk( zhimax(:,:,:), 'T', 1. )
135         
136      ! --- If ice drift field is too fast, use an appropriate time step for advection --- !       
137      zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )    ! CFL test for stability
138      zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
139      IF( lk_mpp )   CALL mpp_max( zcfl )
140     
141      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
142      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
143      ENDIF
144     
145!!      IF( zcfl > 0.5_wp .AND. lwp ) THEN
146!!         ncfl = ncfl + 1
147!!         IF( ncfl > 0 ) THEN   
148!!            WRITE(cltmp,'(i6.1)') ncfl
149!!            CALL ctl_warn( 'ice_adv: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
150!!         ENDIF
151!!      ENDIF
152
153      SELECT CASE ( nn_limadv )
154         
155                       !=============================!
156      CASE ( 0 )       !==  Ultimate-MACHO scheme  ==!                   
157                       !=============================!
158     
159         CALL wrk_alloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
160     
161         IF( kt == nit000 .AND. lwp ) THEN
162            WRITE(numout,*)''
163            WRITE(numout,*)'ice_adv_umx : Ultimate-MACHO advection scheme'
164            WRITE(numout,*)'~~~~~~~~~~~'
165         ENDIF
166         !
167         zdt = rdt_ice / REAL(initad)
168         
169         ! transport
170         zudy(:,:) = u_ice(:,:) * e2u(:,:)
171         zvdx(:,:) = v_ice(:,:) * e1v(:,:)
172         
173         ! define velocity for advection: u*grad(H)
174         DO jj = 2, jpjm1
175            DO ji = fs_2, fs_jpim1
176               IF    ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
177               ELSEIF( u_ice(ji,jj)                  >  0._wp ) THEN   ;   zcu_box(ji,jj) = u_ice(ji-1,jj)
178               ELSE                                                    ;   zcu_box(ji,jj) = u_ice(ji  ,jj)
179               ENDIF
180               
181               IF    ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
182               ELSEIF( v_ice(ji,jj)                  >  0._wp ) THEN   ;   zcv_box(ji,jj) = v_ice(ji,jj-1)
183               ELSE                                                    ;   zcv_box(ji,jj) = v_ice(ji,jj  )
184               ENDIF
185            END DO
186         END DO
187         
188         ! advection
189         DO jt = 1, initad
190            CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) )       ! Open water area
191            DO jl = 1, jpl
192               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) )      ! Ice area
193               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) )      ! Ice  volume
194               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) )    ! Salt content
195               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) )    ! Age content
196               DO jk = 1, nlay_i
197                  CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) )   ! Ice  heat content
198               END DO
199               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume
200               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content
201               ! MV MP 2016
202               IF ( nn_pnd_scheme > 0 ) THEN
203                  CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_ip(:,:,jl) )  ! Melt pond fraction
204                  CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_ip(:,:,jl) )  ! Melt pond volume
205               ENDIF
206               ! END MV MP 2016
207            END DO
208         END DO
209         !
210         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
211         DO jl = 2, jpl
212            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
213         END DO
214         !
215         CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
216         
217                       !=============================!
218      CASE ( -1 )      !==     Prather scheme      ==!                   
219                       !=============================!
220
221         CALL wrk_alloc( jpi,jpj,            zarea )
222         CALL wrk_alloc( jpi,jpj,1,          z0opw )
223         CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
224         CALL wrk_alloc( jpi,jpj,jpl,        z0ap , z0vp )
225         CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
226         
227         IF( kt == nit000 .AND. lwp ) THEN
228            WRITE(numout,*)''
229            WRITE(numout,*)'ice_adv_xy : Prather advection scheme'
230            WRITE(numout,*)'~~~~~~~~~~~'
231         ENDIF
232         
233         zarea(:,:) = e1e2t(:,:)
234         
235         !-------------------------
236         ! transported fields                                       
237         !-------------------------
238         z0opw(:,:,1) = ato_i(:,:) * e1e2t(:,:)             ! Open water area
239         DO jl = 1, jpl
240            z0snw (:,:,jl)  = v_s  (:,:,  jl) * e1e2t(:,:)  ! Snow volume
241            z0ice(:,:,jl)   = v_i  (:,:,  jl) * e1e2t(:,:)  ! Ice  volume
242            z0ai  (:,:,jl)  = a_i  (:,:,  jl) * e1e2t(:,:)  ! Ice area
243            z0smi (:,:,jl)  = smv_i(:,:,  jl) * e1e2t(:,:)  ! Salt content
244            z0oi (:,:,jl)   = oa_i (:,:,  jl) * e1e2t(:,:)  ! Age content
245            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e1e2t(:,:)  ! Snow heat content
246            DO jk = 1, nlay_i
247               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
248            END DO
249            ! MV MP 2016
250            IF ( nn_pnd_scheme > 0 ) THEN
251               z0ap  (:,:,jl)  = a_ip (:,:,jl) * e1e2t(:,:) ! Melt pond fraction
252               z0vp  (:,:,jl)  = v_ip (:,:,jl) * e1e2t(:,:) ! Melt pond volume
253            ENDIF
254            ! END MV MP 2016
255
256         END DO
257
258
259         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
260            DO jt = 1, initad
261               CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
262                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
263               CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
264                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
265               DO jl = 1, jpl
266                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
267                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
268                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
269                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
270                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
271                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
272                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
273                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
274                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
275                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
276                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
277                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
278                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
279                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
280                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
281                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
282                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
283                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
284                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
285                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
286                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
287                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
288                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
289                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
290                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
291                     CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
292                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
293                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
294                     CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
295                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
296                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
297                  END DO
298                  ! MV MP 2016
299                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
300                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
301                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
302                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
303                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   --
304                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
305                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
306                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
307                  ! END MV MP 2016
308               END DO
309            END DO
310         ELSE
311            DO jt = 1, initad
312               CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
313                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
314               CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
315                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
316               DO jl = 1, jpl
317                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
318                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
319                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
320                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
321                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
322                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
323                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
324                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
325                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
326                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
327                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
328                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
329                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
330                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
331                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
332                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
333                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
334                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
335                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
336                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
337                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
338                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
339                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
340                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
341                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
342                     CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
343                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
344                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
345                     CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
346                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
347                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
348                  END DO
349                  ! MV MP 2016
350                  IF ( nn_pnd_scheme > 0 ) THEN
351                     CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
352                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
353                     CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &
354                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
355                     CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   ---
356                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
357                     CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &
358                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
359                  ENDIF
360                  ! END MV MP 2016
361               END DO
362            END DO
363         ENDIF
364
365         !-------------------------------------------
366         ! Recover the properties from their contents
367         !-------------------------------------------
368         ato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:)
369         DO jl = 1, jpl
370            v_i  (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:)
371            v_s  (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:)
372            smv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:)
373            oa_i (:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:)
374            a_i  (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:)
375            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:)
376            DO jk = 1, nlay_i
377               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:)
378            END DO
379            ! MV MP 2016
380            IF ( nn_pnd_scheme > 0 ) THEN
381               a_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:)
382               v_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:)
383            ENDIF
384            ! END MV MP 2016
385         END DO
386
387         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
388         DO jl = 2, jpl
389            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
390         END DO
391         
392         CALL wrk_dealloc( jpi,jpj,            zarea )
393         CALL wrk_dealloc( jpi,jpj,1,          z0opw )
394         CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
395         ! MV MP 2016
396         CALL wrk_dealloc( jpi,jpj,jpl,        z0ap , z0vp )
397         ! END MV MP 2016
398         CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
399
400      END SELECT
401     
402      ! --- diags ---
403      DO jj = 1, jpj
404         DO ji = 1, jpi
405            diag_trp_ei (ji,jj) = ( SUM( e_i  (ji,jj,1:nlay_i,:) ) -  zeiold(ji,jj) ) * r1_rdtice
406            diag_trp_es (ji,jj) = ( SUM( e_s  (ji,jj,1:nlay_s,:) ) -  zesold(ji,jj) ) * r1_rdtice
407            diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:)          ) - zsmvold(ji,jj) ) * r1_rdtice
408            diag_trp_vi (ji,jj) =   SUM(   v_i(ji,jj,:)            -  zviold(ji,jj,:) ) * r1_rdtice
409            diag_trp_vs (ji,jj) =   SUM(   v_s(ji,jj,:)            -  zvsold(ji,jj,:) ) * r1_rdtice
410         END DO
411      END DO
412     
413      IF( nn_limdyn == 2) THEN
414
415         ! zap small areas
416         CALL lim_var_zapsmall
417           
418         !--- Thickness correction in case too high --- !
419         DO jl = 1, jpl
420            DO jj = 1, jpj
421               DO ji = 1, jpi
422                 
423                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
424                     
425                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
426                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
427                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
428                     
429                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 
430                     
431                     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. &
432                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
433                       
434                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
435                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
436                       
437                        ! small correction due to *rswitch for a_i
438                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
439                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
440                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
441                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
442                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
443
444                        ! MV MP 2016
445                        IF ( nn_pnd_scheme > 0 ) THEN
446                           a_ip (ji,jj,jl)        = rswitch * a_ip (ji,jj,jl)
447                           v_ip (ji,jj,jl)        = rswitch * v_ip (ji,jj,jl)
448                        ENDIF
449                        ! END MV MP 2016
450                                               
451                     ENDIF
452                     
453                  ENDIF
454               
455               END DO
456            END DO
457         END DO
458         ! -------------------------------------------------
459         
460         ! Force the upper limit of ht_i to always be < hi_max (99 m).
461         DO jj = 1, jpj
462            DO ji = 1, jpi
463               ! MV MP 2016
464               za_old = a_i(ji,jj,jpl)
465               ! END MV MP 2016
466               rswitch         = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
467               ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
468               a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
469               ! MV MP 2016
470               IF ( nn_pnd_scheme > 0 ) THEN
471                  ! correct pond fraction to avoid a_ip > a_i
472                  a_ip(ji,jj,jpl) = a_ip(ji,jj,jpl) * a_i(ji,jj,jpl) / MAX( za_old , epsi20 ) * rswitch
473               ENDIF
474               ! END MP 2016
475            END DO
476         END DO
477
478
479      ENDIF
480         
481      !------------------------------------------------------------
482      ! Impose a_i < amax if no ridging/rafting or in mono-category
483      !------------------------------------------------------------
484      !
485      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
486      IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2
487         DO jl = 1, jpl
488            DO jj = 1, jpj
489               DO ji = 1, jpi
490                  rswitch       = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) )
491                  zda           = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp )  &
492                     &                    * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 )
493                  a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda
494               END DO
495            END DO
496         END DO
497      ENDIF
498     
499      ! --- agglomerate variables -----------------
500      vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
501      vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
502      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
503
504      ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i)
505      IF ( nn_pnd_scheme > 0 ) THEN
506          at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 )
507          vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 )
508      ENDIF
509      ! END MP 2016
510     
511      ! --- open water = 1 if at_i=0 --------------------------------
512      WHERE( at_i == 0._wp ) ato_i = 1._wp 
513     
514      ! conservation test
515      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
516       
517      ! -------------------------------------------------
518      ! control prints
519      ! -------------------------------------------------
520      IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
521      !
522      IF( nn_timing == 1 )  CALL timing_stop('iceadv')
523      !
524   END SUBROUTINE ice_adv
525
526#else
527   !!----------------------------------------------------------------------
528   !!   Default option         Empty Module                No sea-ice model
529   !!----------------------------------------------------------------------
530CONTAINS
531   SUBROUTINE ice_adv        ! Empty routine
532   END SUBROUTINE ice_adv
533#endif
534   !!======================================================================
535END MODULE iceadv
536
Note: See TracBrowser for help on using the repository browser.