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.
limupdate2.F90 in branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 12.7 KB
Line 
1MODULE limupdate2
2   !!======================================================================
3   !!                     ***  MODULE  limupdate2  ***
4   !!   LIM-3 : Update of sea-ice global variables at the end of the time step
5   !!======================================================================
6   !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
7   !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
8   !!----------------------------------------------------------------------
9#if defined key_lim3
10   !!----------------------------------------------------------------------
11   !!   'key_lim3'                                      LIM3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!    lim_update2   : computes update of sea-ice global variables from trend terms
14   !!----------------------------------------------------------------------
15   USE sbc_oce         ! Surface boundary condition: ocean fields
16   USE sbc_ice         ! Surface boundary condition: ice fields
17   USE dom_ice
18   USE dom_oce
19   USE phycst          ! physical constants
20   USE ice
21   USE thd_ice         ! LIM thermodynamic sea-ice variables
22   USE limitd_th
23   USE limvar
24   USE prtctl          ! Print control
25   USE lbclnk          ! lateral boundary condition - MPP exchanges
26   USE wrk_nemo        ! work arrays
27   USE timing          ! Timing
28   USE limcons         ! conservation tests
29   USE limctl
30   USE lib_mpp         ! MPP library
31   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
32   USE in_out_manager
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   lim_update2   ! routine called by ice_step
38
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE lim_update2( kt )
49      !!-------------------------------------------------------------------
50      !!               ***  ROUTINE lim_update2  ***
51      !!               
52      !! ** Purpose :  Computes update of sea-ice global variables at
53      !!               the end of the time step.
54      !!
55      !!---------------------------------------------------------------------
56      INTEGER, INTENT(in) ::   kt    ! number of iteration
57      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
58      REAL(wp) ::   zsal
59      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
60      !!-------------------------------------------------------------------
61      IF( nn_timing == 1 )  CALL timing_start('limupdate2')
62
63      IF( kt == nit000 .AND. lwp ) THEN
64         WRITE(numout,*) ' lim_update2 '
65         WRITE(numout,*) ' ~~~~~~~~~~~ '
66      ENDIF
67
68      ! conservation test
69      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
70
71      !----------------------------------------------------------------------
72      ! Constrain the thickness of the smallest category above himin
73      !----------------------------------------------------------------------
74      DO jj = 1, jpj 
75         DO ji = 1, jpi
76            rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes
77            ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch
78            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
79               a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin
80               oa_i(ji,jj,1) = oa_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin
81            ENDIF
82         END DO
83      END DO
84     
85      !-----------------------------------------------------
86      ! ice concentration should not exceed amax
87      !-----------------------------------------------------
88      at_i(:,:) = 0._wp
89      DO jl = 1, jpl
90         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
91      END DO
92
93      DO jl  = 1, jpl
94         DO jj = 1, jpj
95            DO ji = 1, jpi
96               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
97                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
98                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
99               ENDIF
100            END DO
101         END DO
102      END DO
103
104      !---------------------
105      ! Ice salinity
106      !---------------------
107      IF (  nn_icesal == 2  ) THEN
108         DO jl = 1, jpl
109            DO jj = 1, jpj 
110               DO ji = 1, jpi
111                  zsal            = smv_i(ji,jj,jl)
112                  ! salinity stays in bounds
113                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
114                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )
115                  ! associated salt flux
116                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
117               END DO
118            END DO
119         END DO
120      ENDIF
121
122      !----------------------------------------------------
123      ! Rebin categories with thickness out of bounds
124      !----------------------------------------------------
125      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )
126
127      !-----------------
128      ! zap small values
129      !-----------------
130      CALL lim_var_zapsmall
131
132      !------------------------------------------------------------------------------
133      ! Corrections to avoid wrong values                                        |
134      !------------------------------------------------------------------------------
135      ! Ice drift
136      !------------
137      DO jj = 2, jpjm1
138         DO ji = 2, jpim1
139            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
140               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
141               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
142               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
143               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
144            ENDIF
145         END DO
146      END DO
147      !lateral boundary conditions
148      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
149      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
150      !mask velocities
151      u_ice(:,:) = u_ice(:,:) * umask(:,:,1)
152      v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)
153 
154      ! -------------------------------------------------
155      ! Diagnostics
156      ! -------------------------------------------------
157      DO jl  = 1, jpl
158         oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice / rday   ! ice natural aging
159         afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
160      END DO
161      afx_tot = afx_thd + afx_dyn
162
163      DO jj = 1, jpj
164         DO ji = 1, jpi           
165            ! heat content variation (W.m-2)
166            diag_heat(ji,jj) = diag_heat(ji,jj) -  &
167               &               ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
168               &                 SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
169               &               ) * r1_rdtice   
170            ! salt, volume
171            diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
172            diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
173            diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
174         END DO
175      END DO
176
177      ! conservation test
178      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
179
180      ! necessary calls (at least for coupling)
181      CALL lim_var_glo2eqv
182      CALL lim_var_agg(2)
183
184      ! -------------------------------------------------
185      ! control prints
186      ! -------------------------------------------------
187      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )   ! control print
188
189      IF(ln_ctl) THEN   ! Control print
190         CALL prt_ctl_info(' ')
191         CALL prt_ctl_info(' - Cell values : ')
192         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
193         CALL prt_ctl(tab2d_1=e1e2t      , clinfo1=' lim_update2  : cell area   :')
194         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :')
195         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :')
196         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update2  : vt_s        :')
197         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update2  : strength    :')
198         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update2  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
199         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update2  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
200
201         DO jl = 1, jpl
202            CALL prt_ctl_info(' ')
203            CALL prt_ctl_info(' - Category : ', ivar1=jl)
204            CALL prt_ctl_info('   ~~~~~~~~~~')
205            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update2  : ht_i        : ')
206            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update2  : ht_s        : ')
207            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update2  : t_su        : ')
208            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : t_snow      : ')
209            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update2  : sm_i        : ')
210            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update2  : o_i         : ')
211            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
212            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : a_i_b       : ')
213            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
214            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_i_b       : ')
215            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
216            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_s_b       : ')
217            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1        : ')
218            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1_b      : ')
219            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2        : ')
220            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2_b      : ')
221            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
222            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow_b    : ')
223            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
224            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update2  : smv_i_b     : ')
225            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
226            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update2  : oa_i_b      : ')
227
228            DO jk = 1, nlay_i
229               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
230               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update2  : t_i       : ')
231            END DO
232         END DO
233
234         CALL prt_ctl_info(' ')
235         CALL prt_ctl_info(' - Heat / FW fluxes : ')
236         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
237         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
238
239         CALL prt_ctl_info(' ')
240         CALL prt_ctl_info(' - Stresses : ')
241         CALL prt_ctl_info('   ~~~~~~~~~~ ')
242         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update2 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
243         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update2 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
244         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update2 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
245      ENDIF
246   
247      IF( nn_timing == 1 )  CALL timing_stop('limupdate2')
248
249   END SUBROUTINE lim_update2
250#else
251   !!----------------------------------------------------------------------
252   !!   Default option         Empty Module               No sea-ice model
253   !!----------------------------------------------------------------------
254CONTAINS
255   SUBROUTINE lim_update2     ! Empty routine
256   END SUBROUTINE lim_update2
257
258#endif
259
260END MODULE limupdate2
Note: See TracBrowser for help on using the repository browser.