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

Last change on this file since 8514 was 8514, checked in by clem, 3 years ago

changes in style - part5 - almost done

File size: 10.3 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   USE iom            !
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.
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, zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc
56      REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag
57      !!----------------------------------------------------------------------
58      IF( nn_timing == 1 )   CALL timing_start('icecor')
59      !
60      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN
61         WRITE(numout,*)
62         WRITE(numout,*) 'ice_cor:  correct sea ice variables if out of bounds ' 
63         WRITE(numout,*) '~~~~~~~'
64      ENDIF
65      !                             !--- conservation test
66      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
67      !
68      !
69      IF( kn == 2 ) THEN
70         !                          !-----------------------------------------------------
71         !                          !  thickness of the smallest category above himin    !
72         !                          !-----------------------------------------------------
73         WHERE( a_i(:,:,1) >= epsi20 )   ;   ht_i(:,:,1) = v_i (:,:,1) / a_i(:,:,1)
74         ELSEWHERE                       ;   ht_i(:,:,1) = 0._wp
75         END WHERE
76         WHERE( ht_i(:,:,1) < rn_himin )     a_i (:,:,1) = a_i (:,:,1) * ht_i(:,:,1) / rn_himin
77         !
78      ENDIF
79      !                             !-----------------------------------------------------
80      !                             !  ice concentration should not exceed amax          !
81      !                             !-----------------------------------------------------
82      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
83      DO jl  = 1, jpl
84         WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:)
85      END DO
86   
87      !                             !-----------------------------------------------------
88      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        !
89      !                             !-----------------------------------------------------
90         zzc = rhoic * r1_rdtice
91         DO jl = 1, jpl
92            DO jj = 1, jpj 
93               DO ji = 1, jpi
94                  zsal = smv_i(ji,jj,jl)
95                  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)  )
96                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux
97               END DO
98            END DO
99         END DO
100      ENDIF
101
102      !                             !-----------------------------------------------------
103      !                             !  Rebin categories with thickness out of bounds     !
104      !                             !-----------------------------------------------------
105      IF ( jpl > 1 )   CALL ice_itd_reb( kt )
106
107      !                             !-----------------------------------------------------
108      CALL ice_var_zapsmall         !  Zap small values                                  !
109      !                             !-----------------------------------------------------
110
111      !                             !-----------------------------------------------------
112      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values !
113         DO jj = 2, jpjm1           !-----------------------------------------------------
114            DO ji = 2, jpim1
115               IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice
116                  IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side
117                  IF ( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side
118                  IF ( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side
119                  IF ( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side
120               ENDIF
121            END DO
122         END DO
123         CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions
124      ENDIF
125
126!!gm I guess the trends are only out on demand
127!!   So please, only do this is it exite an iom_use of on a these variables
128!!   furthermore, only allocate the diag_ arrays in this case
129!!   and do the iom_put here so that it is only a local allocation
130!!gm
131      !                             !-----------------------------------------------------
132      SELECT CASE( kn )             !  Diagnostics                                       !
133      !                             !-----------------------------------------------------
134      CASE( 1 )                        !--- dyn trend diagnostics
135         !
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_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_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         !
155         oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice   ! ice natural aging incrementation
156         !
157!!gm   here I think the number of ice cat is too small to use a SUM instruction...
158         DO jj = 1, jpj
159            DO ji = 1, jpi           
160               !                 ! heat content variation (W.m-2)
161               diag_heat(ji,jj) = diag_heat(ji,jj) - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    & 
162                  &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice
163               !                 ! salt, volume
164               diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
165               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
166               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
167            END DO
168         END DO
169         !                       ! concentration tendency (total + thermo)
170         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice
171         afx_tot(:,:) = afx_tot(:,:) + zafx(:,:)
172         IF( iom_use('afxthd') )   CALL iom_put( 'afxthd' , zafx(:,:) )
173         IF( iom_use('afxtot') )   CALL iom_put( 'afxtot' , afx_tot(:,:) )
174         !
175      END SELECT
176      !
177      !                                !--- conservation test
178      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
179      !
180      !                                !--- control prints
181      IF( ln_ctl )                    CALL ice_prt3D( 'icecor' )
182      IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
183      !
184      IF( nn_timing == 1 )   CALL timing_stop('icecor')
185      !
186   END SUBROUTINE ice_cor
187
188#else
189   !!----------------------------------------------------------------------
190   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
191   !!----------------------------------------------------------------------
192#endif
193
194   !!======================================================================
195END MODULE icecor
Note: See TracBrowser for help on using the repository browser.