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 branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 @ 4649

Last change on this file since 4649 was 4649, checked in by clem, 10 years ago

finalizing LIM3 heat budget conservation + multiple minor bugs corrections
File size: 12.6 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   !!----------------------------------------------------------------------
8#if defined key_lim3
9   !!----------------------------------------------------------------------
10   !!   'key_lim3'                                      LIM3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!    lim_update1   : computes update of sea-ice global variables from trend terms
13   !!----------------------------------------------------------------------
14   USE limrhg          ! ice rheology
15
16   USE dom_oce
17   USE oce             ! dynamics and tracers variables
18   USE in_out_manager
19   USE sbc_oce         ! Surface boundary condition: ocean fields
20   USE sbc_ice         ! Surface boundary condition: ice fields
21   USE dom_ice
22   USE phycst          ! physical constants
23   USE ice
24   USE limdyn
25   USE limtrp
26   USE limthd
27   USE limsbc
28   USE limdiahsb
29   USE limwri
30   USE limrst
31   USE thd_ice         ! LIM thermodynamic sea-ice variables
32   USE par_ice
33   USE limitd_th
34   USE limitd_me
35   USE limvar
36   USE prtctl           ! Print control
37   USE lbclnk           ! lateral boundary condition - MPP exchanges
38   USE wrk_nemo         ! work arrays
39   USE lib_fortran     ! glob_sum
40   USE in_out_manager   ! I/O manager
41   USE iom              ! I/O manager
42   USE lib_mpp          ! MPP library
43   USE timing          ! Timing
44   USE limcons        ! conservation tests
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC   lim_update1   ! routine called by ice_step
50
51      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
52         
53   !! * Substitutions
54#  include "vectopt_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
57   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE lim_update1
63      !!-------------------------------------------------------------------
64      !!               ***  ROUTINE lim_update1  ***
65      !!               
66      !! ** Purpose :  Computes update of sea-ice global variables at
67      !!               the end of the time step.
68      !!               Address pathological cases
69      !!               This place is very important
70      !!               
71      !! ** Method  : 
72      !!    Ice speed from ice dynamics
73      !!    Ice thickness, Snow thickness, Temperatures, Lead fraction
74      !!      from advection and ice thermodynamics
75      !!
76      !! ** Action  : -
77      !!---------------------------------------------------------------------
78      INTEGER ::   ji, jj, jk, jl, jm    ! dummy loop indices
79      INTEGER ::   jbnd1, jbnd2
80      INTEGER ::   i_ice_switch
81      REAL(wp) ::   zinda, zindb, zindsn, zindic
82      REAL(wp) ::   zh, zdvres, zsal
83      REAL(wp) ::   zat_i_old, ztmelts
84      !
85      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
86      !!-------------------------------------------------------------------
87      IF( nn_timing == 1 )  CALL timing_start('limupdate1')
88
89      IF( ln_limdyn ) THEN 
90
91      ! conservation test
92      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
93
94      !-----------------
95      ! zap small values
96      !-----------------
97      CALL lim_itd_me_zapsmall
98
99      CALL lim_var_glo2eqv
100     
101      !----------------------------------------------------
102      ! Rebin categories with thickness out of bounds
103      !----------------------------------------------------
104      DO jm = 1, jpm
105         jbnd1 = ice_cat_bounds(jm,1)
106         jbnd2 = ice_cat_bounds(jm,2)
107         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
108      END DO
109
110      at_i(:,:) = 0._wp
111      DO jl = 1, jpl
112         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
113      END DO
114
115      !----------------------------------------------------
116      ! ice concentration should not exceed amax
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._wp
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      DO jm = 1, jpm
138         jbnd1 = ice_cat_bounds(jm,1)
139         jbnd2 = ice_cat_bounds(jm,2)
140         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
141         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
142         ENDIF
143      END DO
144
145      !-----------------
146      ! zap small values
147      !-----------------
148      CALL lim_itd_me_zapsmall
149
150      !---------------------
151      ! Ice salinity bounds
152      !---------------------
153      IF (  num_sal == 2  ) THEN
154         DO jl = 1, jpl
155            DO jj = 1, jpj 
156               DO ji = 1, jpi
157                  zsal            = smv_i(ji,jj,jl)
158                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
159                  ! salinity stays in bounds
160                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
161                  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) )
162                  ! associated salt flux
163                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
164               END DO
165            END DO
166         END DO
167      ENDIF
168
169      ! -------------------------------------------------
170      ! Diagnostics
171      ! -------------------------------------------------
172      d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:)
173      d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:)
174      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:)
175      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:) 
176      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)   
177      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:) 
178      d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)
179      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:)
180      d_smv_i_trp(:,:,:)   = 0._wp
181      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
182
183      ! conservation test
184      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
185
186      IF(ln_ctl) THEN   ! Control print
187         CALL prt_ctl_info(' ')
188         CALL prt_ctl_info(' - Cell values : ')
189         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
190         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update1  : cell area   :')
191         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :')
192         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :')
193         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update1  : vt_s        :')
194         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update1  : strength    :')
195         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
196         CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update1  : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :')
197         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update1  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
198
199         DO jl = 1, jpl
200            CALL prt_ctl_info(' ')
201            CALL prt_ctl_info(' - Category : ', ivar1=jl)
202            CALL prt_ctl_info('   ~~~~~~~~~~')
203            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update1  : ht_i        : ')
204            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update1  : ht_s        : ')
205            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update1  : t_su        : ')
206            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : t_snow      : ')
207            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update1  : sm_i        : ')
208            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ')
209            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ')
210            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_a_i     : ')
211            CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_a_i_trp   : ')
212            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ')
213            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_i     : ')
214            CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_i_trp   : ')
215            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
216            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_s     : ')
217            CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_s_trp   : ')
218            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ')
219            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i1    : ')
220            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : de_i1_trp   : ')
221            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ')
222            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i2    : ')
223            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : de_i2_trp   : ')
224            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
225            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update1  : old_e_snow  : ')
226            CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : d_e_s_trp   : ')
227            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ')
228            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update1  : old_smv_i   : ')
229            CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update1  : d_smv_i_trp : ')
230            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ')
231            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ')
232            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ')
233
234            DO jk = 1, nlay_i
235               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
236               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update1  : t_i       : ')
237            END DO
238         END DO
239
240         CALL prt_ctl_info(' ')
241         CALL prt_ctl_info(' - Heat / FW fluxes : ')
242         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
243         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
244
245         CALL prt_ctl_info(' ')
246         CALL prt_ctl_info(' - Stresses : ')
247         CALL prt_ctl_info('   ~~~~~~~~~~ ')
248         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update1 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
249         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update1 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
250         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update1 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
251      ENDIF
252   
253      ENDIF ! ln_limdyn
254
255      IF( nn_timing == 1 )  CALL timing_stop('limupdate1')
256   END SUBROUTINE lim_update1
257#else
258   !!----------------------------------------------------------------------
259   !!   Default option         Empty Module               No sea-ice model
260   !!----------------------------------------------------------------------
261CONTAINS
262   SUBROUTINE lim_update1     ! Empty routine
263   END SUBROUTINE lim_update1
264
265#endif
266
267END MODULE limupdate1
Note: See TracBrowser for help on using the repository browser.