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

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90 @ 5055

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

LIM3: removing par_ice and put jpl, nlay_i and nlay_s as namelist parameters
File size: 12.1 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.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_update1   : 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 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   !! * Substitutions
52#  include "vectopt_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
55   !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $
56   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE lim_update1( kt )
61      !!-------------------------------------------------------------------
62      !!               ***  ROUTINE lim_update1  ***
63      !!               
64      !! ** Purpose :  Computes update of sea-ice global variables at
65      !!               the end of the dynamics.
66      !!               
67      !!---------------------------------------------------------------------
68      INTEGER, INTENT(in) :: kt    ! number of iteration
69      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
70      INTEGER  ::   i_ice_switch
71      REAL(wp) ::   zsal
72      !
73      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
74      !!-------------------------------------------------------------------
75      IF( nn_timing == 1 )  CALL timing_start('limupdate1')
76
77      IF( ln_limdyn ) THEN
78
79      IF( kt == nit000 .AND. lwp ) THEN
80         WRITE(numout,*) ' lim_update1 ' 
81         WRITE(numout,*) ' ~~~~~~~~~~~ '
82      ENDIF
83
84      ! conservation test
85      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
86
87      !-----------------
88      ! zap small values
89      !-----------------
90      CALL lim_var_zapsmall
91
92      CALL lim_var_glo2eqv
93     
94      !----------------------------------------------------
95      ! Rebin categories with thickness out of bounds
96      !----------------------------------------------------
97      IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl)
98
99      at_i(:,:) = 0._wp
100      DO jl = 1, jpl
101         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
102      END DO
103
104      !----------------------------------------------------
105      ! ice concentration should not exceed amax
106      !-----------------------------------------------------
107      DO jl  = 1, jpl
108         DO jj = 1, jpj
109            DO ji = 1, jpi
110               IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
111                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) )
112                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
113               ENDIF
114            END DO
115         END DO
116      END DO
117
118      at_i(:,:) = 0._wp
119      DO jl = 1, jpl
120         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
121      END DO
122   
123      ! --------------------------------------
124      ! Final thickness distribution rebinning
125      ! --------------------------------------
126      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl)
127
128      !-----------------
129      ! zap small values
130      !-----------------
131      CALL lim_var_zapsmall
132
133      !---------------------
134      ! Ice salinity bounds
135      !---------------------
136      IF (  num_sal == 2  ) THEN
137         DO jl = 1, jpl
138            DO jj = 1, jpj 
139               DO ji = 1, jpi
140                  zsal            = smv_i(ji,jj,jl)
141                  smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl)
142                  ! salinity stays in bounds
143                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
144                  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) )
145                  ! associated salt flux
146                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
147               END DO
148            END DO
149         END DO
150      ENDIF
151
152      ! -------------------------------------------------
153      ! Diagnostics
154      ! -------------------------------------------------
155      DO jl  = 1, jpl
156         afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
157      END DO
158
159      d_u_ice_dyn(:,:)     = u_ice(:,:)     - u_ice_b(:,:)
160      d_v_ice_dyn(:,:)     = v_ice(:,:)     - v_ice_b(:,:)
161      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - a_i_b  (:,:,:)
162      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - v_s_b  (:,:,:) 
163      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - v_i_b  (:,:,:)   
164      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - e_s_b  (:,:,:,:) 
165      d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:)
166      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - oa_i_b (:,:,:)
167      d_smv_i_trp(:,:,:)   = 0._wp
168      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:)
169
170      ! conservation test
171      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
172
173      IF(ln_ctl) THEN   ! Control print
174         CALL prt_ctl_info(' ')
175         CALL prt_ctl_info(' - Cell values : ')
176         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
177         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update1  : cell area   :')
178         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :')
179         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :')
180         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update1  : vt_s        :')
181         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update1  : strength    :')
182         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
183         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 :')
184         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update1  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
185
186         DO jl = 1, jpl
187            CALL prt_ctl_info(' ')
188            CALL prt_ctl_info(' - Category : ', ivar1=jl)
189            CALL prt_ctl_info('   ~~~~~~~~~~')
190            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update1  : ht_i        : ')
191            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update1  : ht_s        : ')
192            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update1  : t_su        : ')
193            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : t_snow      : ')
194            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update1  : sm_i        : ')
195            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ')
196            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ')
197            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : a_i_b       : ')
198            CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_a_i_trp   : ')
199            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ')
200            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_i_b       : ')
201            CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_i_trp   : ')
202            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
203            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ')
204            CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_s_trp   : ')
205            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ')
206            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1_b      : ')
207            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : de_i1_trp   : ')
208            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ')
209            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2_b      : ')
210            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : de_i2_trp   : ')
211            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
212            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ')
213            CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : d_e_s_trp   : ')
214            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ')
215            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update1  : smv_i_b     : ')
216            CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update1  : d_smv_i_trp : ')
217            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ')
218            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update1  : oa_i_b      : ')
219            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ')
220
221            DO jk = 1, nlay_i
222               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
223               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update1  : t_i       : ')
224            END DO
225         END DO
226
227         CALL prt_ctl_info(' ')
228         CALL prt_ctl_info(' - Heat / FW fluxes : ')
229         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
230         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
231
232         CALL prt_ctl_info(' ')
233         CALL prt_ctl_info(' - Stresses : ')
234         CALL prt_ctl_info('   ~~~~~~~~~~ ')
235         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update1 : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
236         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update1 : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
237         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update1 : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
238      ENDIF
239   
240      ENDIF ! ln_limdyn
241
242      IF( nn_timing == 1 )  CALL timing_stop('limupdate1')
243   END SUBROUTINE lim_update1
244#else
245   !!----------------------------------------------------------------------
246   !!   Default option         Empty Module               No sea-ice model
247   !!----------------------------------------------------------------------
248CONTAINS
249   SUBROUTINE lim_update1     ! Empty routine
250   END SUBROUTINE lim_update1
251
252#endif
253
254END MODULE limupdate1
Note: See TracBrowser for help on using the repository browser.