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 @ 8491

Last change on this file since 8491 was 8491, checked in by clem, 6 years ago

debug icevar.F90
File size: 11.8 KB
Line 
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        ! ocean domain
16   USE phycst         ! physical constants
17   USE ice            ! sea-ice: variable
18   USE ice1D          ! sea-ice: thermodynamic sea-ice variables
19   USE iceitd         ! sea-ice: rebining
20   USE icevar         ! sea-ice: operations
21   USE icectl         ! sea-ice: 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   ! called by icestp.F90
33
34   !! * Substitutions
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
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 corrections on sea-ice global variables at
48      !!              the end of the dynamics.
49      !!----------------------------------------------------------------------
50      INTEGER, INTENT(in) ::   kt    ! number of iteration
51      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo
52      !
53      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
54      REAL(wp) ::   zsal, zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc
55      !!----------------------------------------------------------------------
56      IF( nn_timing == 1 )   CALL timing_start('icecor')
57      !
58      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN
59         WRITE(numout,*)
60         WRITE(numout,*)' icecor :  correct sea ice variables if out of bounds ' 
61         WRITE(numout,*)' ~~~~~~~'
62      ENDIF
63      !                             !--- conservation test
64      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
65      !
66      !
67      !                             !-----------------------------------------------------
68      IF( kn == 2 ) THEN            !  thickness of the smallest category above himin    !
69         !                          !-----------------------------------------------------
70         !
71         DO jj = 1, jpj 
72            DO ji = 1, jpi
73!!gm replace this
74               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes
75               ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch
76!!gm by more readable coding (not slower coding since already a IF in the loop):
77!               IF( a_i(ji,jj,1) >= epsi20 )   ht_i(ji,jj,1) = v_i (ji,jj,1) / a_i(ji,jj,1)
78!!gm
79               IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
80                  a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin
81               ENDIF
82            END DO
83         END DO
84         !
85      ENDIF
86
87      !                             !-----------------------------------------------------
88      at_i(:,:) = a_i(:,:,1)        !  ice concentration should not exceed amax          !
89      DO jl = 2, jpl                !-----------------------------------------------------
90         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
91      END DO
92      !
93!!gm   Question   it seams to me that we have the following equality (dropping the "(ji,jj)":
94!      ( 1. - ( 1. - rn_amax_2d / at_i ) ) =  ( 1. - ( at_i - rn_amax_2d ) / at_i )
95!                                          =  ( at_i - ( at_i - rn_amax_2d ) ) / at_i
96!                                          =  ( + rn_amax_2d  ) / at_i
97!                                          =  rn_amax_2d / at_i
98!     No ?  if yes see "!!gm   better" juste below
99!gm
100      DO jl  = 1, jpl
101         DO jj = 1, jpj
102            DO ji = 1, jpi
103               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
104                  a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
105!!gm   better:    a_i(ji,jj,jl) = a_i(ji,jj,jl) * rn_amax_2d(ji,jj) / at_i(ji,jj)           
106               ENDIF
107            END DO
108         END DO
109      END DO
110!!gm  Other question:  why testing a_i(ji,jj,jl) > 0._wp ?   a_i is >=0, a multiplication by 0 does not change the results....
111!!gm  so at the end, the loop can be recoded without IF as:
112!      WHERE( at_i(:,:) > rn_amax_2d(:,:) )
113!         DO jl  = 1, jpl
114!            a_i(:,:,jl) = a_i(:,:,jl) * MAX( rn_amax_2d(:,:), at_i(:,:) ) / at_i(:,:)
115!         END DO
116!      END WHERE
117!!gm  No?
118   
119      !                             !-----------------------------------------------------
120      IF (  nn_icesal == 2  ) THEN  !  Ice salinity bounds                               !
121      !                             !-----------------------------------------------------
122         zzc = rhoic * r1_rdtice
123         DO jl = 1, jpl                  ! salinity stays in bounds [Simin,Simax]
124            DO jj = 1, jpj 
125               DO ji = 1, jpi
126                  IF( v_i(ji,jj,jl) > 0._wp ) THEN   ! clem: useless IF ???
127                     zsal = smv_i(ji,jj,jl)
128                     smv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , smv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  )
129                     ! associated salt flux
130                     sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc
131                  ENDIF
132               END DO
133            END DO
134         END DO
135      ENDIF
136
137      !                             !-----------------------------------------------------
138      !                             !  Rebin categories with thickness out of bounds     !
139      !                             !-----------------------------------------------------
140      IF ( jpl > 1 )   CALL ice_itd_reb
141
142      !                             !-----------------------------------------------------
143      CALL ice_var_zapsmall         !  Zap small values                                  !
144      !                             !-----------------------------------------------------
145
146      !                             !-----------------------------------------------------
147      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values !
148         DO jj = 2, jpjm1           !-----------------------------------------------------
149            DO ji = 2, jpim1
150               IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice
151                  IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side
152                  IF ( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side
153                  IF ( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side
154                  IF ( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side
155               ENDIF
156            END DO
157         END DO
158         CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions
159         !
160!!gm  I think masking here is unnecessary, u_ice already masked and we only introduce zeros in the field
161         u_ice(:,:) = u_ice(:,:) * umask(:,:,1)                            ! mask velocities
162         v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)
163      ENDIF
164
165!!gm I guess the trends are only out on demand
166!!   So please, only do this is it exite an iom_use of on a these variables
167!!   furthermore, only allocate the diag_ arrays in this case
168!!   and do the iom_put here so that it is only a local allocation
169!!gm
170      !                             !-----------------------------------------------------
171      SELECT CASE( kn )             !  Diagnostics                                       !
172      !                             !-----------------------------------------------------
173      CASE( 1 )                        !--- dyn trend diagnostics
174         DO jl  = 1, jpl
175            afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
176         END DO
177         !
178!!gm   here I think the number of ice cat is too small to use a SUM instruction...
179         DO jj = 1, jpj
180            DO ji = 1, jpi           
181               !                 ! heat content variation (W.m-2)
182               diag_heat(ji,jj) = - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    & 
183                  &                  + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice
184               !                 ! salt, volume
185               diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
186               diag_vice(ji,jj) = SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
187               diag_vsnw(ji,jj) = SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
188            END DO
189         END DO
190         !
191      CASE( 2 )                        !--- thermo trend diagnostics & ice aging
192         !
193         DO jl  = 1, jpl
194            oa_i(:,:,jl) = oa_i(:,:,jl) +   a_i(:,:,jl)                   * rdt_ice       ! ice natural aging incrementation
195            afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice     ! thermo tendancy
196         END DO
197         afx_tot(:,:) = afx_thd(:,:) + afx_dyn(:,:)
198         !
199!!gm   here I think the number of ice cat is too small to use a SUM instruction...
200         DO jj = 1, jpj
201            DO ji = 1, jpi           
202               !                 ! heat content variation (W.m-2)
203               diag_heat(ji,jj) = diag_heat(ji,jj) - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    & 
204                  &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice
205               !                 ! salt, volume
206               diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
207               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
208               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
209            END DO
210         END DO
211         !
212      END SELECT
213      !
214      !                                !--- conservation test
215      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
216      !
217      !                                !--- control prints
218      IF( ln_ctl )                    CALL ice_prt3D( 'icecor' )
219      IF( ln_limctl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
220      !
221      IF( nn_timing == 1 )   CALL timing_stop('icecor')
222      !
223   END SUBROUTINE ice_cor
224
225#else
226   !!----------------------------------------------------------------------
227   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
228   !!----------------------------------------------------------------------
229#endif
230
231   !!======================================================================
232END MODULE icecor
Note: See TracBrowser for help on using the repository browser.