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

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90

Last change on this file was 8892, checked in by frrh, 7 years ago

Commit updates with debugging write statements.
File size: 10.6 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( nn_timing == 1 )   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      IF( kn == 2 ) THEN
69         !                          !-----------------------------------------------------
70         !                          !  thickness of the smallest category above himin    !
71         !                          !-----------------------------------------------------
72         WHERE( a_i(:,:,1) >= epsi20 )   ;   h_i(:,:,1) = v_i (:,:,1) / a_i(:,:,1)
73         ELSEWHERE                       ;   h_i(:,:,1) = 0._wp
74         END WHERE
75         WHERE( h_i(:,:,1) < rn_himin )      a_i(:,:,1) = a_i (:,:,1) * h_i(:,:,1) / rn_himin
76         !
77      ENDIF
78      !                             !-----------------------------------------------------
79      !                             !  ice concentration should not exceed amax          !
80      !                             !-----------------------------------------------------
81      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
82      DO jl  = 1, jpl
83         WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:)
84      END DO
85   
86      !                             !-----------------------------------------------------
87      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        !
88      !                             !-----------------------------------------------------
89         zzc = rhoic * r1_rdtice
90         DO jl = 1, jpl
91            DO jj = 1, jpj 
92               DO ji = 1, jpi
93                  zsal = sv_i(ji,jj,jl)
94                  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)  )
95                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux
96               END DO
97            END DO
98         END DO
99      ENDIF
100
101      !                             !-----------------------------------------------------
102      !                             !  Rebin categories with thickness out of bounds     !
103      !                             !-----------------------------------------------------
104      IF ( jpl > 1 )   CALL ice_itd_reb( kt )
105write(numout,*) "RSRH out of ice_itd_reb in icecor", kt ; flush(numout)
106      !                             !-----------------------------------------------------
107      CALL ice_var_zapsmall         !  Zap small values                                  !
108      !                             !-----------------------------------------------------
109
110      !                             !-----------------------------------------------------
111      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values !
112         DO jj = 2, jpjm1           !-----------------------------------------------------
113            DO ji = 2, jpim1
114               IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice
115                  IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side
116                  IF ( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side
117                  IF ( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side
118                  IF ( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side
119               ENDIF
120            END DO
121         END DO
122         CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions
123      ENDIF
124
125!!gm I guess the trends are only out on demand
126!!   So please, only do this is it exite an iom_use of on a these variables
127!!   furthermore, only allocate the diag_ arrays in this case
128!!   and do the iom_put here so that it is only a local allocation
129!!gm
130      !                             !-----------------------------------------------------
131      SELECT CASE( kn )             !  Diagnostics                                       !
132      !                             !-----------------------------------------------------
133      CASE( 1 )                        !--- dyn trend diagnostics
134         !
135write(numout,*) "RSRH GG icecor CASE 1", kt ; flush(numout)
136!!gm   here I think the number of ice cat is too small to use a SUM instruction...
137         DO jj = 1, jpj
138            DO ji = 1, jpi           
139               !                 ! heat content variation (W.m-2)
140               diag_heat(ji,jj) = - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    & 
141                  &                  + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice
142               !                 ! salt, volume
143               diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
144               diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice
145               diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice
146            END DO
147         END DO
148         !                       ! concentration tendency (dynamics)
149         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
150         afx_tot(:,:) = zafx(:,:)
151         IF( iom_use('afxdyn') )   CALL iom_put( 'afxdyn' , zafx(:,:) )
152         !
153      CASE( 2 )                        !--- thermo trend diagnostics & ice aging
154         !
155write(numout,*) "RSRH GG icecor CASE 2", kt ; flush(numout)
156         oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice   ! ice natural aging incrementation
157         !
158!!gm   here I think the number of ice cat is too small to use a SUM instruction...
159         DO jj = 1, jpj
160            DO ji = 1, jpi           
161               !                 ! heat content variation (W.m-2)
162               diag_heat(ji,jj) = diag_heat(ji,jj) - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    & 
163                  &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice
164               !                 ! salt, volume
165               diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
166               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice
167               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice
168            END DO
169         END DO
170         !                       ! concentration tendency (total + thermo)
171         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice
172         afx_tot(:,:) = afx_tot(:,:) + zafx(:,:)
173         IF( iom_use('afxthd') )   CALL iom_put( 'afxthd' , zafx(:,:) )
174         IF( iom_use('afxtot') )   CALL iom_put( 'afxtot' , afx_tot(:,:) )
175         !
176      END SELECT
177write(numout,*) "RSRH HH icecor", kt ; flush(numout)
178      !
179      ! controls
180      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
181      IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints
182      IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                   ! prints
183      IF( nn_timing == 1 )   CALL timing_stop ('icecor')                                                             ! timing
184write(numout,*) "RSRH leaving icecor", kt ; flush(numout)
185      !
186   END SUBROUTINE ice_cor
187
188#else
189   !!----------------------------------------------------------------------
190   !!   Default option           Dummy module         NO ESIM sea-ice model
191   !!----------------------------------------------------------------------
192#endif
193
194   !!======================================================================
195END MODULE icecor
Note: See TracBrowser for help on using the repository browser.