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.
icecor.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90 @ 8426

Last change on this file since 8426 was 8426, checked in by clem, 7 years ago

last routine names to be changed
File size: 8.7 KB
RevLine 
[8426]1MODULE icecor
2   !!======================================================================
3   !!                     ***  MODULE  icecor  ***
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   !!    ice_cor   : computes update of sea-ice global variables from trend terms
14   !!----------------------------------------------------------------------
15   USE dom_oce
16   USE phycst          ! physical constants
17   USE ice
18   USE ice1D           ! LIM thermodynamic sea-ice variables
19   USE iceitd
20   USE icevar
21   USE icectl          ! control prints
22   !
23   USE in_out_manager  ! I/O manager
24   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
25   USE lbclnk          ! lateral boundary condition - MPP link
26   USE lib_mpp         ! MPP library
27   USE timing          ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   ice_cor
33
34   !! * Substitutions
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
38   !! $Id: icecor.F90 8378 2017-07-26 13:55:59Z clem $
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE ice_cor( kt , kn )
44      !!-------------------------------------------------------------------
45      !!               ***  ROUTINE ice_cor  ***
46      !!               
47      !! ** Purpose :  Computes update of sea-ice global variables at
48      !!               the end of the dynamics.
49      !!               
50      !!---------------------------------------------------------------------
51      INTEGER, INTENT(in) ::   kt    ! number of iteration
52      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo
53      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
54      REAL(wp) ::   zsal
55      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
56      !!-------------------------------------------------------------------
57      IF( nn_timing == 1 )  CALL timing_start('icecor')
58
59      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN
60         WRITE(numout,*)
61         WRITE(numout,*)' icecor ' 
62         WRITE(numout,*)' ~~~~~~ '
63      ENDIF
64
65      ! conservation test
66      IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
67
68      !----------------------------------------------------------------------
69      ! Constrain the thickness of the smallest category above himin
70      !----------------------------------------------------------------------
71      IF( kn == 2 ) THEN
72         
73         DO jj = 1, jpj 
74            DO ji = 1, jpi
75               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes
76               ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch
77               IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
78                  a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin
79               ENDIF
80            END DO
81         END DO
82
83      ENDIF
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_2d(ji,jj) .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_2d(ji,jj) / at_i(ji,jj) ) )
98               ENDIF
99            END DO
100         END DO
101      END DO
102   
103      !---------------------
104      ! Ice salinity bounds
105      !---------------------
106      IF (  nn_icesal == 2  ) THEN
107         DO jl = 1, jpl
108            DO jj = 1, jpj 
109               DO ji = 1, jpi
110                  zsal            = smv_i(ji,jj,jl)
111                  ! salinity stays in bounds
112                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
113                  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) )
114                  ! associated salt flux
115                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
116               END DO
117            END DO
118         END DO
119      ENDIF
120
121      !----------------------------------------------------
122      ! Rebin categories with thickness out of bounds
123      !----------------------------------------------------
124      IF ( jpl > 1 ) CALL ice_itd_reb
125
126      !-----------------
127      ! zap small values
128      !-----------------
129      CALL ice_var_zapsmall
130
131      !----------------------------------------------
132      ! Ice drift. Corrections to avoid wrong values
133      !----------------------------------------------
134      IF( kn == 2 ) THEN
135         DO jj = 2, jpjm1
136            DO ji = 2, jpim1
137               IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
138                  IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
139                  IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
140                  IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
141                  IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
142               ENDIF
143            END DO
144         END DO
145         !lateral boundary conditions
146         CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )
147         !mask velocities
148         u_ice(:,:) = u_ice(:,:) * umask(:,:,1)
149         v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)
150      ENDIF
151     
152      ! -------------------------------------------------
153      ! Diagnostics
154      ! -------------------------------------------------
155      IF( kn == 1 ) THEN
156         DO jl  = 1, jpl
157            afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
158         END DO
159         
160         DO jj = 1, jpj
161            DO ji = 1, jpi           
162               ! heat content variation (W.m-2)
163               diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
164                  &                   SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
165                  &                 ) * r1_rdtice
166               ! salt, volume
167               diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
168               diag_vice(ji,jj) = SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
169               diag_vsnw(ji,jj) = SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
170            END DO
171         END DO
172         
173      ELSEIF( kn == 2 ) THEN
174         
175         DO jl  = 1, jpl
176            oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice          ! ice natural aging
177            afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
178         END DO
179         afx_tot = afx_thd + afx_dyn
180         
181         DO jj = 1, jpj
182            DO ji = 1, jpi           
183               ! heat content variation (W.m-2)
184               diag_heat(ji,jj) = diag_heat(ji,jj) -  &
185                  &               ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
186                  &                 SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
187                  &               ) * r1_rdtice   
188               ! salt, volume
189               diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
190               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
191               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
192            END DO
193         END DO
194         
195      ENDIF
196     
197      ! conservation test
198      IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
199
200      ! control prints
201      IF( ln_ctl )       CALL ice_prt3D( 'icecor' )
202      IF( ln_limctl .AND. kn == 2 )  &
203         &               CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
204   
205      IF( nn_timing == 1 )  CALL timing_stop('icecor')
206
207   END SUBROUTINE ice_cor
208
209#endif
210
211END MODULE icecor
Note: See TracBrowser for help on using the repository browser.