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

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 5053

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

LIM3 cleaning continues. No change in the physics besides the introduction of the monocategory sea ice
File size: 14.1 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.6  !  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 limrhg          ! ice rheology
16
17   USE dom_oce
18   USE oce             ! dynamics and tracers variables
19   USE in_out_manager
20   USE sbc_oce         ! Surface boundary condition: ocean fields
21   USE sbc_ice         ! Surface boundary condition: ice fields
22   USE dom_ice
23   USE phycst          ! physical constants
24   USE ice
25   USE limdyn
26   USE limtrp
27   USE limthd
28   USE limsbc
29   USE limdiahsb
30   USE limwri
31   USE limrst
32   USE thd_ice         ! LIM thermodynamic sea-ice variables
33   USE par_ice
34   USE limitd_th
35   USE limitd_me
36   USE limvar
37   USE prtctl           ! Print control
38   USE lbclnk           ! lateral boundary condition - MPP exchanges
39   USE wrk_nemo         ! work arrays
40   USE lib_fortran     ! glob_sum
41   USE timing          ! Timing
42   USE limcons        ! conservation tests
43   USE limctl
44
45   IMPLICIT NONE
46   PRIVATE
47
48   PUBLIC   lim_update2   ! routine called by ice_step
49
50   !! * Substitutions
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
54   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
55   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE lim_update2( kt )
60      !!-------------------------------------------------------------------
61      !!               ***  ROUTINE lim_update2  ***
62      !!               
63      !! ** Purpose :  Computes update of sea-ice global variables at
64      !!               the end of the time step.
65      !!
66      !!---------------------------------------------------------------------
67      INTEGER, INTENT(in) :: kt    ! number of iteration
68      INTEGER  ::   ji, jj, jk, jl    ! dummy loop indices
69      INTEGER  ::   i_ice_switch
70      REAL(wp) ::   zh, zsal
71      !
72      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
73      !!-------------------------------------------------------------------
74      IF( nn_timing == 1 )  CALL timing_start('limupdate2')
75
76      IF( kt == nit000 .AND. lwp ) THEN
77         WRITE(numout,*) ' lim_update2 '
78         WRITE(numout,*) ' ~~~~~~~~~~~ '
79      ENDIF
80
81      ! conservation test
82      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
83
84      !-----------------
85      ! zap small values
86      !-----------------
87      CALL lim_var_agg( 1 )
88      CALL lim_var_zapsmall
89      CALL lim_var_glo2eqv
90
91      !----------------------------------------------------
92      ! Rebin categories with thickness out of bounds
93      !----------------------------------------------------
94      IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl)
95
96      !----------------------------------------------------------------------
97      ! Constrain the thickness of the smallest category above hiclim
98      !----------------------------------------------------------------------
99      DO jj = 1, jpj 
100         DO ji = 1, jpi
101            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < hiclim ) THEN
102               zh             = hiclim / ht_i(ji,jj,1)
103               ht_s(ji,jj,1) = ht_s(ji,jj,1) * zh
104               ht_i(ji,jj,1) = ht_i(ji,jj,1) * zh
105               a_i (ji,jj,1) = a_i(ji,jj,1)  / zh
106            ENDIF
107         END DO
108      END DO
109     
110      !-----------------------------------------------------
111      ! ice concentration should not exceed amax
112      !-----------------------------------------------------
113      at_i(:,:) = 0._wp
114      DO jl = 1, jpl
115         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
116      END DO
117
118      DO jl  = 1, jpl
119         DO jj = 1, jpj
120            DO ji = 1, jpi
121               IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
122                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) )
123                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
124               ENDIF
125            END DO
126         END DO
127      END DO
128
129      at_i(:,:) = 0.0
130      DO jl = 1, jpl
131         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
132      END DO
133
134      ! --------------------------------------
135      ! Final thickness distribution rebinning
136      ! --------------------------------------
137      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )
138
139      !-----------------
140      ! zap small values
141      !-----------------
142      CALL lim_var_zapsmall
143
144      !---------------------
145      ! 2.11) Ice salinity
146      !---------------------
147      IF (  num_sal == 2  ) THEN
148         DO jl = 1, jpl
149            DO jj = 1, jpj 
150               DO ji = 1, jpi
151                  zsal            = smv_i(ji,jj,jl)
152                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
153                  ! salinity stays in bounds
154                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
155                  smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
156                  ! associated salt flux
157                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
158               END DO ! ji
159            END DO ! jj
160         END DO !jl
161      ENDIF
162
163      !------------------------------------------------------------------------------
164      ! 2) Corrections to avoid wrong values                                        |
165      !------------------------------------------------------------------------------
166      ! Ice drift
167      !------------
168      DO jj = 2, jpjm1
169         DO ji = 2, jpim1
170            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
171               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
172               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
173               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
174               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
175            ENDIF
176         END DO
177      END DO
178      !lateral boundary conditions
179      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
180      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
181      !mask velocities
182      u_ice(:,:) = u_ice(:,:) * tmu(:,:)
183      v_ice(:,:) = v_ice(:,:) * tmv(:,:)
184 
185      ! for outputs
186      CALL lim_var_glo2eqv            ! equivalent variables (outputs)
187      CALL lim_var_agg(2)             ! aggregate ice thickness categories
188
189      ! -------------------------------------------------
190      ! Diagnostics
191      ! -------------------------------------------------
192      DO jl  = 1, jpl
193         afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
194      END DO
195      afx_tot = afx_thd + afx_dyn
196
197      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - a_i_b(:,:,:) 
198      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - v_s_b(:,:,:)
199      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - v_i_b(:,:,:) 
200      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - e_s_b(:,:,:,:) 
201      d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:)
202      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - oa_i_b (:,:,:)
203      d_smv_i_thd(:,:,:) = 0._wp
204      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:)
205
206      ! heat content variation (W.m-2)
207      DO jj = 1, jpj
208         DO ji = 1, jpi           
209            diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) +  & 
210               &                     SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) )    &
211               &                   ) * unit_fac * r1_rdtice / area(ji,jj)   
212         END DO
213      END DO
214
215      ! conservation test
216      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
217
218      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print
219
220      IF(ln_ctl) THEN   ! Control print
221         CALL prt_ctl_info(' ')
222         CALL prt_ctl_info(' - Cell values : ')
223         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
224         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update2  : cell area   :')
225         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :')
226         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :')
227         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update2  : vt_s        :')
228         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update2  : strength    :')
229         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update2  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
230         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update2  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
231
232         DO jl = 1, jpl
233            CALL prt_ctl_info(' ')
234            CALL prt_ctl_info(' - Category : ', ivar1=jl)
235            CALL prt_ctl_info('   ~~~~~~~~~~')
236            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update2  : ht_i        : ')
237            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update2  : ht_s        : ')
238            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update2  : t_su        : ')
239            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : t_snow      : ')
240            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update2  : sm_i        : ')
241            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update2  : o_i         : ')
242            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
243            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : a_i_b       : ')
244            CALL prt_ctl(tab2d_1=d_a_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_a_i_thd   : ')
245            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
246            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_i_b       : ')
247            CALL prt_ctl(tab2d_1=d_v_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_i_thd   : ')
248            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
249            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_s_b       : ')
250            CALL prt_ctl(tab2d_1=d_v_s_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_s_thd   : ')
251            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1        : ')
252            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1_b      : ')
253            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : de_i1_thd   : ')
254            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2        : ')
255            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2_b      : ')
256            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : de_i2_thd   : ')
257            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
258            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow_b    : ')
259            CALL prt_ctl(tab2d_1=d_e_s_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : d_e_s_thd   : ')
260            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
261            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update2  : smv_i_b     : ')
262            CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl)        , clinfo1= ' lim_update2  : d_smv_i_thd : ')
263            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
264            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update2  : oa_i_b      : ')
265            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ')
266
267            DO jk = 1, nlay_i
268               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
269               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update2  : t_i       : ')
270            END DO
271         END DO
272
273         CALL prt_ctl_info(' ')
274         CALL prt_ctl_info(' - Heat / FW fluxes : ')
275         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
276         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
277
278         CALL prt_ctl_info(' ')
279         CALL prt_ctl_info(' - Stresses : ')
280         CALL prt_ctl_info('   ~~~~~~~~~~ ')
281         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update2 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
282         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update2 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
283         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update2 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
284      ENDIF
285   
286      IF( nn_timing == 1 )  CALL timing_stop('limupdate2')
287
288   END SUBROUTINE lim_update2
289#else
290   !!----------------------------------------------------------------------
291   !!   Default option         Empty Module               No sea-ice model
292   !!----------------------------------------------------------------------
293CONTAINS
294   SUBROUTINE lim_update2     ! Empty routine
295   END SUBROUTINE lim_update2
296
297#endif
298
299END MODULE limupdate2
Note: See TracBrowser for help on using the repository browser.