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

source: branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 5350

Last change on this file since 5350 was 5350, checked in by hadcv, 9 years ago

Update to head of the trunk (r5344).

  • Property svn:keywords set to Id
File size: 26.1 KB
RevLine 
[825]1MODULE limtrp
2   !!======================================================================
3   !!                       ***  MODULE limtrp   ***
4   !! LIM transport ice model : sea-ice advection/diffusion
5   !!======================================================================
[2715]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   !!----------------------------------------------------------------------
[825]10#if defined key_lim3
11   !!----------------------------------------------------------------------
[834]12   !!   'key_lim3'                                      LIM3 sea-ice model
[825]13   !!----------------------------------------------------------------------
14   !!   lim_trp      : advection/diffusion process of sea ice
15   !!----------------------------------------------------------------------
[3625]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
[5123]23   USE limvar         !
24   !
[3625]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) 
[5123]31   USE timing         ! Timing
[4688]32   USE limcons        ! conservation tests
[5123]33   USE limctl         ! control prints
[825]34
35   IMPLICIT NONE
36   PRIVATE
37
[5123]38   PUBLIC   lim_trp    ! called by sbcice_lim
[825]39
[5123]40   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2 
41
[825]42   !! * Substitution
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
[4161]45   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]46   !! $Id$
[2715]47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[825]48   !!----------------------------------------------------------------------
49CONTAINS
50
[921]51   SUBROUTINE lim_trp( kt ) 
[825]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      !!---------------------------------------------------------------------
[5123]63      INTEGER, INTENT(in) ::   kt           ! number of iteration
[2715]64      !
[5123]65      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
[2715]66      INTEGER  ::   initad                  ! number of sub-timestep for the advection
[4990]67      REAL(wp) ::   zcfl , zusnit           !   -      -
[5123]68      CHARACTER(len=80) ::   cltmp
[2715]69      !
[5134]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
[5123]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) ::    zdv, zvi, zvs, zsmv, zes, zei
78      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[2715]79      !!---------------------------------------------------------------------
[4161]80      IF( nn_timing == 1 )  CALL timing_start('limtrp')
[825]81
[5350]82      CALL wrk_alloc( jpi,jpj,            zsm, zatold, zeiold, zesold )
83      CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
84      CALL wrk_alloc( jpi,jpj,1,          z0opw )
85      CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
86      CALL wrk_alloc( jpi,jpj,jpl,        zhimax, zviold, zvsold, zsmvold )
[825]87
[2715]88      IF( numit == nstart .AND. lwp ) THEN
89         WRITE(numout,*)
90         IF( ln_limdyn ) THEN   ;   WRITE(numout,*) 'lim_trp : Ice transport '
91         ELSE                   ;   WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
92         ENDIF
93         WRITE(numout,*) '~~~~~~~~~~~~'
[5123]94         ncfl = 0                ! nb of time step with CFL > 1/2
[2715]95      ENDIF
[5123]96
97      zsm(:,:) = e12t(:,:)
[2715]98     
99      !                             !-------------------------------------!
100      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
101         !                          !-------------------------------------!
[4688]102
103         ! conservation test
[5123]104         IF( ln_limdiahsb )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4688]105
[5123]106         ! mass and salt flux init
[4161]107         zviold(:,:,:)  = v_i(:,:,:)
[5123]108         zvsold(:,:,:)  = v_s(:,:,:)
109         zsmvold(:,:,:) = smv_i(:,:,:)
110         zeiold(:,:)    = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 
111         zesold(:,:)    = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 
[825]112
[5123]113         !--- Thickness correction init. -------------------------------
114         zatold(:,:) = SUM( a_i(:,:,:), dim=3 )
[5350]115         DO jl = 1, jpl
116            DO jj = 1, jpj
117               DO ji = 1, jpi
118                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
119                  ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
120                  ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
121               END DO
122            END DO
123         END DO
[4161]124         !---------------------------------------------------------------------
[5350]125         ! Record max of the surrounding ice thicknesses for correction
[4161]126         ! in case advection creates ice too thick.
127         !---------------------------------------------------------------------
[5134]128         zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
[4161]129         DO jl = 1, jpl
130            DO jj = 2, jpjm1
131               DO ji = 2, jpim1
[5134]132                  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) )
[4161]133               END DO
134            END DO
135            CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
136         END DO
137         
[5123]138         !=============================!
139         !==      Prather scheme     ==!
140         !=============================!
141
142         ! If ice drift field is too fast, use an appropriate time step for advection.         
143         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )         ! CFL test for stability
144         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
145         IF(lk_mpp )   CALL mpp_max( zcfl )
146
147         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
148         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
149         ENDIF
150
151         IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1
[5350]152!!         IF( lwp ) THEN
153!!            IF( ncfl > 0 ) THEN   
154!!               WRITE(cltmp,'(i6.1)') ncfl
155!!               CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
156!!            ELSE
157!!            !  WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 '
158!!            ENDIF
159!!         ENDIF
[5123]160
[825]161         !-------------------------
[2715]162         ! transported fields                                       
[825]163         !-------------------------
[5134]164         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area
[2715]165         DO jl = 1, jpl
[5134]166            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
167            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
168            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
169            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
170            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
171            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
[5123]172            DO jk = 1, nlay_i
[5134]173               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
[5123]174            END DO
[825]175         END DO
176
[921]177
[2715]178         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
[5123]179            DO jt = 1, initad
[5134]180               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
[5123]181                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
[5134]182               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
[5123]183                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
[825]184               DO jl = 1, jpl
[5134]185                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
[2715]186                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
[5134]187                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
[2715]188                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
[5134]189                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
[2715]190                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
[5134]191                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
[2715]192                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
[5134]193                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
[2715]194                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
[5134]195                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
[2715]196                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
[5134]197                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
[2715]198                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
[5134]199                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
[2715]200                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
[5134]201                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
[2715]202                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
[5134]203                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
[2715]204                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
[5134]205                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
[2715]206                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
[5134]207                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
[2715]208                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
[5123]209                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
[5134]210                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
[4870]211                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
212                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[5134]213                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
[4870]214                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
215                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]216                  END DO
217               END DO
218            END DO
219         ELSE
[5123]220            DO jt = 1, initad
[5134]221               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
[5123]222                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
[5134]223               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
[5123]224                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
[825]225               DO jl = 1, jpl
[5134]226                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
[2715]227                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
[5134]228                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
[2715]229                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
[5134]230                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
[2715]231                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
[5134]232                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
[2715]233                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
[5134]234                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
[2715]235                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
[5134]236                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
[2715]237                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
[5134]238                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
[2715]239                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
[5134]240                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
[2715]241                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
[5134]242                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
[2715]243                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
[5134]244                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &
[2715]245                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
[5134]246                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
[2715]247                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
[5134]248                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
[2715]249                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
[4870]250                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
[5134]251                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
[4870]252                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
253                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[5134]254                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
[4870]255                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
256                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]257                  END DO
258               END DO
259            END DO
260         ENDIF
261
262         !-------------------------------------------
263         ! Recover the properties from their contents
264         !-------------------------------------------
[5134]265         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
[825]266         DO jl = 1, jpl
[5134]267            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
268            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
269            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
270            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
271            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
272            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
[5123]273            DO jk = 1, nlay_i
[5134]274               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
[5123]275            END DO
[825]276         END DO
277
[5123]278         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
279         DO jl = 2, jpl
280            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
281         END DO
282
[921]283         !------------------------------------------------------------------------------!
[5123]284         ! Diffusion of Ice fields                 
[921]285         !------------------------------------------------------------------------------!
[825]286
[5123]287         !
[2715]288         !--------------------------------
289         !  diffusion of open water area
290         !--------------------------------
291         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
292         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
293            DO ji = 1 , fs_jpim1   ! vector opt.
[5123]294               pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
295                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
296               pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
297                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
[2715]298            END DO
299         END DO
300         !
[5134]301         CALL lim_hdf( ato_i (:,:) )
[2715]302
[921]303         !------------------------------------
[2715]304         !  Diffusion of other ice variables
[921]305         !------------------------------------
[825]306         DO jl = 1, jpl
[2715]307         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
308            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
309               DO ji = 1 , fs_jpim1   ! vector opt.
[5123]310                  pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   &
311                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
312                  pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   &
313                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
[825]314               END DO
[921]315            END DO
[825]316
[5123]317            CALL lim_hdf( v_i  (:,:,  jl) )
318            CALL lim_hdf( v_s  (:,:,  jl) )
319            CALL lim_hdf( smv_i(:,:,  jl) )
320            CALL lim_hdf( oa_i (:,:,  jl) )
321            CALL lim_hdf( a_i  (:,:,  jl) )
322            CALL lim_hdf( e_s  (:,:,1,jl) )
[825]323            DO jk = 1, nlay_i
[5123]324               CALL lim_hdf( e_i(:,:,jk,jl) )
[2715]325            END DO
326         END DO
[825]327
[921]328         !------------------------------------------------------------------------------!
[5123]329         ! limit ice properties after transport                           
[921]330         !------------------------------------------------------------------------------!
[5123]331!!gm & cr   :  MAX should not be active if adv scheme is positive !
[825]332         DO jl = 1, jpl
333            DO jj = 1, jpj
334               DO ji = 1, jpi
[5123]335                  v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) )
336                  v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) )
337                  smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) )
338                  oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) )
339                  a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) )
340                  e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) )
[825]341               END DO
342            END DO
343
[4688]344            DO jk = 1, nlay_i
345               DO jj = 1, jpj
346                  DO ji = 1, jpi
[5123]347                     e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) )
348                  END DO
349               END DO
350            END DO
351         END DO
352!!gm & cr
[4688]353
[5350]354         ! --- diags ---
355         DO jj = 1, jpj
356            DO ji = 1, jpi
357               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice
358               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice
359
360               diag_trp_vi (ji,jj) = SUM(   v_i(ji,jj,:) -  zviold(ji,jj,:) ) * r1_rdtice
361               diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice
362               diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice
363            END DO
364         END DO
365
[5123]366         ! zap small areas
367         CALL lim_var_zapsmall
368
369         !--- Thickness correction in case too high --------------------------------------------------------
[4161]370         DO jl = 1, jpl
371            DO jj = 1, jpj
372               DO ji = 1, jpi
373
374                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
[5134]375
[5350]376                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
377                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
378                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
379                     
[4688]380                     zvi  = v_i  (ji,jj,jl)
381                     zvs  = v_s  (ji,jj,jl)
382                     zsmv = smv_i(ji,jj,jl)
383                     zes  = e_s  (ji,jj,1,jl)
[4990]384                     zei  = SUM( e_i(ji,jj,1:nlay_i,jl) )
[5134]385
386                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 
387
388                     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. &
[5350]389                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
[5134]390
391                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
392                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
393
394                        ! small correction due to *rswitch for a_i
395                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
396                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
397                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
398                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
399                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
400
401                        ! Update mass fluxes
402                        wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
403                        wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
404                        sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
405                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
406                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
407
[4161]408                     ENDIF
409
410                  ENDIF
[5123]411
[825]412               END DO
413            END DO
414         END DO
[4688]415         ! -------------------------------------------------
[5123]416         
417         !--------------------------------------
418         ! Impose a_i < amax in mono-category
419         !--------------------------------------
420         !
421         IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
422            DO jj = 1, jpj
423               DO ji = 1, jpi
424                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax )
425               END DO
426            END DO
427         ENDIF
[825]428
[4990]429         ! --- agglomerate variables -----------------
[4688]430         vt_i (:,:) = 0._wp
431         vt_s (:,:) = 0._wp
432         at_i (:,:) = 0._wp
[825]433         DO jl = 1, jpl
434            DO jj = 1, jpj
435               DO ji = 1, jpi
[5134]436                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl)
437                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl)
438                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
[4688]439               END DO
440            END DO
441         END DO
[825]442
[5134]443         ! --- open water = 1 if at_i=0 --------------------------------
[4161]444         DO jj = 1, jpj
445            DO ji = 1, jpi
[4990]446               rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
[5123]447               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
[4161]448            END DO
[4688]449         END DO     
[4161]450
[4688]451         ! conservation test
452         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]453
[825]454      ENDIF
455
[5123]456      ! -------------------------------------------------
457      ! control prints
458      ! -------------------------------------------------
[5128]459      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
[2715]460      !
[5350]461      CALL wrk_dealloc( jpi,jpj,            zsm, zatold, zeiold, zesold )
462      CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
463      CALL wrk_dealloc( jpi,jpj,1,          z0opw )
464      CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
465      CALL wrk_dealloc( jpi,jpj,jpl,        zviold, zvsold, zhimax, zsmvold )
[5123]466      !
[4161]467      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
[5123]468
[825]469   END SUBROUTINE lim_trp
470
471#else
472   !!----------------------------------------------------------------------
473   !!   Default option         Empty Module                No sea-ice model
474   !!----------------------------------------------------------------------
475CONTAINS
476   SUBROUTINE lim_trp        ! Empty routine
477   END SUBROUTINE lim_trp
478#endif
479   !!======================================================================
480END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.