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.
limupdate1.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 @ 5134

Last change on this file since 5134 was 5134, checked in by clem, 9 years ago

LIM3: changes to constrain ice thickness after advection. Ice volume and concentration are advected while ice thickness is deduced. This could result in overly high thicknesses associated with very low concentrations (in the 5th category)
File size: 10.3 KB
Line 
1MODULE limupdate1
2   !!======================================================================
3   !!                     ***  MODULE  limupdate1  ***
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_update1   : 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 wrk_nemo        ! work arrays
26   USE timing          ! Timing
27   USE limcons         ! conservation tests
28   USE lib_mpp         ! MPP library
29   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30   USE in_out_manager  ! I/O manager
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   lim_update1
36
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
41   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE lim_update1( kt )
47      !!-------------------------------------------------------------------
48      !!               ***  ROUTINE lim_update1  ***
49      !!               
50      !! ** Purpose :  Computes update of sea-ice global variables at
51      !!               the end of the dynamics.
52      !!               
53      !!---------------------------------------------------------------------
54      INTEGER, INTENT(in) ::   kt    ! number of iteration
55      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
56      REAL(wp) ::   zsal
57      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
58      !!-------------------------------------------------------------------
59      IF( nn_timing == 1 )  CALL timing_start('limupdate1')
60
61      IF( ln_limdyn ) THEN
62
63      IF( kt == nit000 .AND. lwp ) THEN
64         WRITE(numout,*) ' lim_update1 ' 
65         WRITE(numout,*) ' ~~~~~~~~~~~ '
66      ENDIF
67
68      ! conservation test
69      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
70
71      CALL lim_var_glo2eqv
72      !----------------------------------------------------
73      ! ice concentration should not exceed amax
74      !-----------------------------------------------------
75      at_i(:,:) = 0._wp
76      DO jl = 1, jpl
77         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
78      END DO
79
80      DO jl  = 1, jpl
81         DO jj = 1, jpj
82            DO ji = 1, jpi
83               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
84                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
85               ENDIF
86            END DO
87         END DO
88      END DO
89   
90      !----------------------------------------------------
91      ! Rebin categories with thickness out of bounds
92      !----------------------------------------------------
93      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl)
94
95      !-----------------
96      ! zap small values
97      !-----------------
98      CALL lim_var_zapsmall
99
100      !---------------------
101      ! Ice salinity bounds
102      !---------------------
103      IF (  nn_icesal == 2  ) THEN
104         DO jl = 1, jpl
105            DO jj = 1, jpj 
106               DO ji = 1, jpi
107                  zsal            = smv_i(ji,jj,jl)
108                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
109                  ! salinity stays in bounds
110                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
111                  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) )
112                  ! associated salt flux
113                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
114               END DO
115            END DO
116         END DO
117      ENDIF
118
119      ! conservation test
120      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
121
122      ! -------------------------------------------------
123      ! Diagnostics
124      ! -------------------------------------------------
125      DO jl  = 1, jpl
126         afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
127      END DO
128
129      ! heat content variation (W.m-2)
130      DO jj = 1, jpj
131         DO ji = 1, jpi           
132            diag_heat_dhc(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
133               &                       SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
134               &                     ) * r1_rdtice   
135         END DO
136      END DO
137
138      ! -------------------------------------------------
139      ! control prints
140      ! -------------------------------------------------
141      IF(ln_ctl) THEN   ! Control print
142         CALL prt_ctl_info(' ')
143         CALL prt_ctl_info(' - Cell values : ')
144         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
145         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update1  : cell area   :')
146         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :')
147         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :')
148         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update1  : vt_s        :')
149         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update1  : strength    :')
150         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
151         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update1  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
152
153         DO jl = 1, jpl
154            CALL prt_ctl_info(' ')
155            CALL prt_ctl_info(' - Category : ', ivar1=jl)
156            CALL prt_ctl_info('   ~~~~~~~~~~')
157            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update1  : ht_i        : ')
158            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update1  : ht_s        : ')
159            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update1  : t_su        : ')
160            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : t_snow      : ')
161            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update1  : sm_i        : ')
162            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ')
163            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ')
164            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : a_i_b       : ')
165            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ')
166            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_i_b       : ')
167            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
168            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ')
169            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1        : ')
170            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1_b      : ')
171            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2        : ')
172            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2_b      : ')
173            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
174            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ')
175            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ')
176            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update1  : smv_i_b     : ')
177            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ')
178            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update1  : oa_i_b      : ')
179
180            DO jk = 1, nlay_i
181               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
182               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update1  : t_i       : ')
183            END DO
184         END DO
185
186         CALL prt_ctl_info(' ')
187         CALL prt_ctl_info(' - Heat / FW fluxes : ')
188         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
189         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
190
191         CALL prt_ctl_info(' ')
192         CALL prt_ctl_info(' - Stresses : ')
193         CALL prt_ctl_info('   ~~~~~~~~~~ ')
194         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update1 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
195         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update1 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
196         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update1 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
197      ENDIF
198   
199      ENDIF ! ln_limdyn
200
201      IF( nn_timing == 1 )  CALL timing_stop('limupdate1')
202   END SUBROUTINE lim_update1
203#else
204   !!----------------------------------------------------------------------
205   !!   Default option         Empty Module               No sea-ice model
206   !!----------------------------------------------------------------------
207CONTAINS
208   SUBROUTINE lim_update1     ! Empty routine
209   END SUBROUTINE lim_update1
210
211#endif
212
213END MODULE limupdate1
Note: See TracBrowser for help on using the repository browser.