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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90 @ 9415

Last change on this file since 9415 was 9415, checked in by clem, 5 years ago

debug ice diffusion: avoid too small ice thicknesses by including a threshold (rn_himin) after ice dynamics. This should prevent the ice-ocean conduction flux from being abnormally strong

File size: 10.3 KB
Line 
1MODULE icecor
2   !!======================================================================
3   !!                     ***  MODULE  icecor  ***
4   !!   ESIM : Corrections on sea-ice 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'                                       ESIM sea-ice model
12   !!----------------------------------------------------------------------
13   !!    ice_cor      : corrections on sea-ice variables
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 iom            ! I/O manager library
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
27   USE lbclnk         ! lateral boundary conditions (or mpp links)
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_cor   ! called by icestp.F90
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
39   !! $Id: icecor.F90 8378 2017-07-26 13:55:59Z clem $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE ice_cor( kt, kn )
45      !!----------------------------------------------------------------------
46      !!               ***  ROUTINE ice_cor  ***
47      !!               
48      !! ** Purpose :   Computes corrections on sea-ice global variables at
49      !!              the end of the dynamics (kn=1) and thermodynamics (kn=2)
50      !!----------------------------------------------------------------------
51      INTEGER, INTENT(in) ::   kt    ! number of iteration
52      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo
53      !
54      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
55      REAL(wp) ::   zsal, zzc
56      REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag
57      !!----------------------------------------------------------------------
58      ! controls
59      IF( ln_timing    )   CALL timing_start('icecor')                                                             ! timing
60      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
61      !
62      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN
63         WRITE(numout,*)
64         WRITE(numout,*) 'ice_cor:  correct sea ice variables if out of bounds ' 
65         WRITE(numout,*) '~~~~~~~'
66      ENDIF
67      !
68      !                             !-----------------------------------------------------
69      !                             !  ice thickness must exceed himin (for ice diff)    !
70      !                             !-----------------------------------------------------
71      WHERE( a_i(:,:,:) >= epsi20 )   ;   h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
72      ELSEWHERE                       ;   h_i(:,:,:) = 0._wp
73      END WHERE
74      WHERE( h_i(:,:,:) < rn_himin )      a_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:) / rn_himin
75      !
76      !                             !-----------------------------------------------------
77      !                             !  ice concentration should not exceed amax          !
78      !                             !-----------------------------------------------------
79      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
80      DO jl  = 1, jpl
81         WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:)
82      END DO
83   
84      !                             !-----------------------------------------------------
85      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        !
86      !                             !-----------------------------------------------------
87         zzc = rhoic * r1_rdtice
88         DO jl = 1, jpl
89            DO jj = 1, jpj 
90               DO ji = 1, jpi
91                  zsal = sv_i(ji,jj,jl)
92                  sv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  )
93                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux
94               END DO
95            END DO
96         END DO
97      ENDIF
98
99      !                             !-----------------------------------------------------
100      !                             !  Rebin categories with thickness out of bounds     !
101      !                             !-----------------------------------------------------
102      IF ( jpl > 1 )   CALL ice_itd_reb( kt )
103
104      !                             !-----------------------------------------------------
105      CALL ice_var_zapsmall         !  Zap small values                                  !
106      !                             !-----------------------------------------------------
107
108      !                             !-----------------------------------------------------
109      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values !
110         DO jj = 2, jpjm1           !-----------------------------------------------------
111            DO ji = 2, jpim1
112               IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice
113                  IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side
114                  IF ( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side
115                  IF ( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side
116                  IF ( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side
117               ENDIF
118            END DO
119         END DO
120         CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions
121      ENDIF
122
123!!gm I guess the trends are only out on demand
124!!   So please, only do this is it exite an iom_use of on a these variables
125!!   furthermore, only allocate the diag_ arrays in this case
126!!   and do the iom_put here so that it is only a local allocation
127!!gm
128      !                             !-----------------------------------------------------
129      SELECT CASE( kn )             !  Diagnostics                                       !
130      !                             !-----------------------------------------------------
131      CASE( 1 )                        !--- dyn trend diagnostics
132         !
133!!gm   here I think the number of ice cat is too small to use a SUM instruction...
134         DO jj = 1, jpj
135            DO ji = 1, jpi           
136               !                 ! heat content variation (W.m-2)
137               diag_heat(ji,jj) = - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    & 
138                  &                  + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice
139               !                 ! salt, volume
140               diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
141               diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice
142               diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice
143            END DO
144         END DO
145         !                       ! concentration tendency (dynamics)
146         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
147         afx_tot(:,:) = zafx(:,:)
148         IF( iom_use('afxdyn') )   CALL iom_put( 'afxdyn' , zafx(:,:) )
149         !
150      CASE( 2 )                        !--- thermo trend diagnostics & ice aging
151         !
152         oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice   ! ice natural aging incrementation
153         !
154!!gm   here I think the number of ice cat is too small to use a SUM instruction...
155         DO jj = 1, jpj
156            DO ji = 1, jpi           
157               !                 ! heat content variation (W.m-2)
158               diag_heat(ji,jj) = diag_heat(ji,jj) - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    & 
159                  &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice
160               !                 ! salt, volume
161               diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
162               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice
163               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice
164            END DO
165         END DO
166         !                       ! concentration tendency (total + thermo)
167         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice
168         afx_tot(:,:) = afx_tot(:,:) + zafx(:,:)
169         IF( iom_use('afxthd') )   CALL iom_put( 'afxthd' , zafx(:,:) )
170         IF( iom_use('afxtot') )   CALL iom_put( 'afxtot' , afx_tot(:,:) )
171         !
172      END SELECT
173      !
174      ! controls
175      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
176      IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints
177      IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                   ! prints
178      IF( ln_timing      )   CALL timing_stop ('icecor')                                                             ! timing
179      !
180   END SUBROUTINE ice_cor
181
182#else
183   !!----------------------------------------------------------------------
184   !!   Default option           Dummy module         NO ESIM sea-ice model
185   !!----------------------------------------------------------------------
186#endif
187
188   !!======================================================================
189END MODULE icecor
Note: See TracBrowser for help on using the repository browser.