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

source: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 7506

Last change on this file since 7506 was 7506, checked in by vancop, 8 years ago

Commit a first set of modifications

  • Property svn:keywords set to Id
File size: 32.9 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 dom_ice        ! ice domain
20   USE ice            ! ice variables
21   USE limadv         ! ice advection
22   USE limhdf         ! ice horizontal diffusion
23   USE limvar         !
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, a second order Prather scheme is used. 
60      !!
61      !! ** action :
62      !!---------------------------------------------------------------------
63      INTEGER, INTENT(in) ::   kt           ! number of iteration
64      !
65      INTEGER  ::   ji, jj, jk, jm , jl, jt      ! dummy loop indices
66      INTEGER  ::   initad                  ! number of sub-timestep for the advection
67      REAL(wp) ::   zcfl , zusnit           !   -      -
68      CHARACTER(len=80) ::   cltmp
69      !
70      REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm
71      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
72      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
73      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
74      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold, zsmvold  ! old ice volume...
75      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness
76      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies
77      REAL(wp), POINTER, DIMENSION(:,:,:)             ::   zhdfptab
78      REAL(wp) ::    zdv, zvi, zvs, zsmv, zes, zei
79      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
80      !!---------------------------------------------------------------------
81      INTEGER                                ::  ihdf_vars  = 6  !!Number of variables in which we apply horizontal diffusion
82                                                                   !!  inside limtrp for each ice category , not counting the
83                                                                   !!  variables corresponding to ice_layers
84      !!---------------------------------------------------------------------
85      IF( nn_timing == 1 )  CALL timing_start('limtrp')
86
87      CALL wrk_alloc( jpi,jpj,            zsm, zatold, zeiold, zesold )
88      CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
89      CALL wrk_alloc( jpi,jpj,1,          z0opw )
90      CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
91      CALL wrk_alloc( jpi,jpj,jpl,        zhimax, zviold, zvsold, zsmvold )
92      CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1,zhdfptab)
93
94      IF( numit == nstart .AND. lwp ) THEN
95         WRITE(numout,*)
96         IF( ln_limdyn ) THEN   ;   WRITE(numout,*) 'lim_trp : Ice transport '
97         ELSE                   ;   WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
98         ENDIF
99         WRITE(numout,*) '~~~~~~~~~~~~'
100         ncfl = 0                ! nb of time step with CFL > 1/2
101      ENDIF
102
103      zsm(:,:) = e12t(:,:)
104     
105      !                             !-------------------------------------!
106      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
107         !                          !-------------------------------------!
108
109         ! conservation test
110         IF( ln_limdiahsb )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
111
112         ! mass and salt flux init
113         zviold(:,:,:)  = v_i(:,:,:)
114         zvsold(:,:,:)  = v_s(:,:,:)
115         zsmvold(:,:,:) = smv_i(:,:,:)
116         zeiold(:,:)    = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 
117         zesold(:,:)    = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 
118
119         !--- Thickness correction init. -------------------------------
120         zatold(:,:) = SUM( a_i(:,:,:), dim=3 )
121         DO jl = 1, jpl
122            DO jj = 1, jpj
123               DO ji = 1, jpi
124                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
125                  ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
126                  ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
127               END DO
128            END DO
129         END DO
130         !---------------------------------------------------------------------
131         ! Record max of the surrounding ice thicknesses for correction
132         ! in case advection creates ice too thick.
133         !---------------------------------------------------------------------
134         zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
135         DO jl = 1, jpl
136            DO jj = 2, jpjm1
137               DO ji = 2, jpim1
138                  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) )
139               END DO
140            END DO
141            CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
142         END DO
143         
144         !=============================!
145         !==      Prather scheme     ==!
146         !=============================!
147
148         ! If ice drift field is too fast, use an appropriate time step for advection.         
149         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )         ! CFL test for stability
150         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
151         IF(lk_mpp )   CALL mpp_max( zcfl )
152
153         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
154         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
155         ENDIF
156
157         IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1
158!!         IF( lwp ) THEN
159!!            IF( ncfl > 0 ) THEN   
160!!               WRITE(cltmp,'(i6.1)') ncfl
161!!               CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
162!!            ELSE
163!!            !  WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 '
164!!            ENDIF
165!!         ENDIF
166
167         !-------------------------
168         ! transported fields                                       
169         !-------------------------
170         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area
171         DO jl = 1, jpl
172            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
173            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
174            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
175            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
176            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
177            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
178           DO jk = 1, nlay_i
179               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
180            END DO
181         END DO
182
183
184         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
185            DO jt = 1, initad
186               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
187                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
188               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
189                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
190               DO jl = 1, jpl
191                  ! SIMIP mass transport diags
192                  diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
193                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
194                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
195                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
196                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
197                  diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
198
199                  diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
200                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume
201                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
202                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume
203                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
204                  diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
205                  ! END SIMIP mass transport diags
206
207!                 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
208!                    &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
209!                 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
210!                    &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
211!                 CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
212!                    &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
213!                 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
214!                    &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
215                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
216                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
217                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
218                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
219                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
220                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
221                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
222                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
223                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
224                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
225                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
226                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
227                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
228                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
229                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
230                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
231                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
232                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
233                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
234                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
235                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
236                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
237                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
238                  END DO
239               END DO
240            END DO
241         ELSE
242            DO jt = 1, initad
243               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
244                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
245               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
246                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
247               DO jl = 1, jpl
248                  diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
249                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
250                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
251                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
252                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
253                  diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
254
255                  diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
256                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
257                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
258                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
259                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
260                  diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
261
262!                 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
263!                    &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
264!                 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
265!                    &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
266!                 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
267!                    &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
268!                 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
269!                    &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
270
271                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
272                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
273                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
274                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
275                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
276                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
277                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
278                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
279                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
280                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
281                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &
282                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
283                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
284                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
285                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
286                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
287                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
288                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
289                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
290                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
291                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
292                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
293                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
294                  END DO
295               END DO
296            END DO
297         ENDIF
298
299         ! SIMIP diags
300         diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) / ( rdt_ice * zusnit )
301         diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) / ( rdt_ice * zusnit )
302
303         !-------------------------------------------
304         ! Recover the properties from their contents
305         !-------------------------------------------
306         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
307         DO jl = 1, jpl
308            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
309            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
310            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
311            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
312            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
313            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
314            DO jk = 1, nlay_i
315               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
316            END DO
317         END DO
318
319         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
320         DO jl = 2, jpl
321            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
322         END DO
323
324         !------------------------------------------------------------------------------!
325         ! Diffusion of Ice fields                 
326         !------------------------------------------------------------------------------!
327         !------------------------------------
328         !  Diffusion of other ice variables
329         !------------------------------------
330         jm=1
331         DO jl = 1, jpl
332         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
333         !   DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
334         !      DO ji = 1 , fs_jpim1   ! vector opt.
335         !         pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   &
336         !            &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
337         !         pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   &
338         !            &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
339         !      END DO
340         !   END DO
341            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
342               DO ji = 1 , fs_jpim1   ! vector opt.
343                  pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,  jl ) ) ) )   &
344                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj)
345                  pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,  jj,  jl ) ) ) )   &
346                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,  jj+1,jl ) ) ) ) * ahiv(ji,jj)
347               END DO
348            END DO
349
350            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1   
351            zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1
352            zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1 
353            zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1
354            zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
355            zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
356         ! Sample of adding more variables to apply lim_hdf using lim_hdf optimization---
357         !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1 
358         !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1
359         !
360         ! and in this example the parameter ihdf_vars musb be changed to 8 (necessary for allocation)
361         !----------------------------------------------------------------------------------------
362            DO jk = 1, nlay_i
363              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
364            END DO
365         END DO
366         !
367         !--------------------------------
368         !  diffusion of open water area
369         !--------------------------------
370         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
371         !DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
372         !   DO ji = 1 , fs_jpim1   ! vector opt.
373         !      pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
374         !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
375         !      pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
376         !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
377         !   END DO
378         !END DO
379         
380         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
381            DO ji = 1 , fs_jpim1   ! vector opt.
382               pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
383                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
384               pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
385                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
386            END DO
387         END DO
388         !
389         zhdfptab(:,:,jm)= ato_i  (:,:);
390         CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i) 
391
392         jm=1
393         DO jl = 1, jpl
394            a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1     
395            v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
396            v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
397            smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
398            oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
399            e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1 
400         ! Sample of adding more variables to apply lim_hdf---------
401         !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
402         !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
403         !-----------------------------------------------------------
404            DO jk = 1, nlay_i
405               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 
406            END DO
407         END DO
408
409         ato_i  (:,:) = zhdfptab(:,:,jm)
410
411         !------------------------------------------------------------------------------!
412         ! limit ice properties after transport                           
413         !------------------------------------------------------------------------------!
414!!gm & cr   :  MAX should not be active if adv scheme is positive !
415         DO jl = 1, jpl
416            DO jj = 1, jpj
417               DO ji = 1, jpi
418                  v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) )
419                  v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) )
420                  smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) )
421                  oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) )
422                  a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) )
423                  e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) )
424               END DO
425            END DO
426
427            DO jk = 1, nlay_i
428               DO jj = 1, jpj
429                  DO ji = 1, jpi
430                     e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) )
431                  END DO
432               END DO
433            END DO
434         END DO
435!!gm & cr
436
437         ! --- diags ---
438         DO jj = 1, jpj
439            DO ji = 1, jpi
440               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice
441               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice
442
443               diag_trp_vi (ji,jj) = SUM(   v_i(ji,jj,:) -  zviold(ji,jj,:) ) * r1_rdtice
444               diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice
445               diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice
446
447               ! SIMIP diagnostics
448               diag_dms_dyn(ji,jj) = rhosn * diag_trp_vs(ji,jj)
449               diag_dmi_dyn(ji,jj) = rhoic * diag_trp_vi(ji,jj)
450            END DO
451         END DO
452
453         ! zap small areas
454         CALL lim_var_zapsmall
455
456         !--- Thickness correction in case too high --------------------------------------------------------
457         DO jl = 1, jpl
458            DO jj = 1, jpj
459               DO ji = 1, jpi
460
461                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
462
463                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
464                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
465                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
466                     
467                     zvi  = v_i  (ji,jj,jl)
468                     zvs  = v_s  (ji,jj,jl)
469                     zsmv = smv_i(ji,jj,jl)
470                     zes  = e_s  (ji,jj,1,jl)
471                     zei  = SUM( e_i(ji,jj,1:nlay_i,jl) )
472
473                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 
474
475                     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. &
476                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
477
478                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
479                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
480
481                        ! small correction due to *rswitch for a_i
482                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
483                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
484                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
485                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
486                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
487
488                        ! Update mass fluxes
489                        wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
490                        wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
491                        sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
492                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
493                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
494
495                        ! SIMIP diagnostic
496                        diag_dms_dyn(ji,jj) = diag_dms_dyn(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
497
498                     ENDIF
499
500                  ENDIF
501
502               END DO
503            END DO
504         END DO
505         ! -------------------------------------------------
506         
507         !--------------------------------------
508         ! Impose a_i < amax in mono-category
509         !--------------------------------------
510         !
511         IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
512            DO jj = 1, jpj
513               DO ji = 1, jpi
514                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax_2d(ji,jj) )
515               END DO
516            END DO
517         ENDIF
518
519         ! --- agglomerate variables -----------------
520         vt_i (:,:) = 0._wp
521         vt_s (:,:) = 0._wp
522         at_i (:,:) = 0._wp
523         DO jl = 1, jpl
524            DO jj = 1, jpj
525               DO ji = 1, jpi
526                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl)
527                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl)
528                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
529               END DO
530            END DO
531         END DO
532
533         ! --- open water = 1 if at_i=0 --------------------------------
534         DO jj = 1, jpj
535            DO ji = 1, jpi
536               rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
537               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
538            END DO
539         END DO     
540
541         ! conservation test
542         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
543
544      ENDIF
545
546      ! -------------------------------------------------
547      ! control prints
548      ! -------------------------------------------------
549      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
550      !
551      CALL wrk_dealloc( jpi,jpj,            zsm, zatold, zeiold, zesold )
552      CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
553      CALL wrk_dealloc( jpi,jpj,1,          z0opw )
554      CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
555      CALL wrk_dealloc( jpi,jpj,jpl,        zviold, zvsold, zhimax, zsmvold )
556      CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars+nlay_i)+1,zhdfptab)
557      !
558      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
559
560   END SUBROUTINE lim_trp
561
562#else
563   !!----------------------------------------------------------------------
564   !!   Default option         Empty Module                No sea-ice model
565   !!----------------------------------------------------------------------
566CONTAINS
567   SUBROUTINE lim_trp        ! Empty routine
568   END SUBROUTINE lim_trp
569#endif
570   !!======================================================================
571END MODULE limtrp
572
Note: See TracBrowser for help on using the repository browser.