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.
limupdate2.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 9.1 KB
Line 
1MODULE limupdate2
2   !!======================================================================
3   !!                     ***  MODULE  limupdate2  ***
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   !!    lim_update2   : computes update of sea-ice global variables from trend terms
14   !!----------------------------------------------------------------------
15   USE sbc_oce         ! Surface boundary condition: ocean fields
16   USE sbc_ice         ! Surface boundary condition: ice fields
17   USE dom_oce
18   USE phycst          ! physical constants
19   USE ice
20   USE thd_ice         ! LIM thermodynamic sea-ice variables
21   USE limitd_th
22   USE limvar
23   USE lbclnk          ! lateral boundary condition - MPP exchanges
24   USE wrk_nemo        ! work arrays
25   USE timing          ! Timing
26   USE limcons         ! conservation tests
27   USE limctl
28   USE lib_mpp         ! MPP library
29   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30   USE in_out_manager
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   lim_update2   ! routine called by ice_step
36
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE lim_update2( kt )
47      !!-------------------------------------------------------------------
48      !!               ***  ROUTINE lim_update2  ***
49      !!               
50      !! ** Purpose :  Computes update of sea-ice global variables at
51      !!               the end of the time step.
52      !!
53      !!---------------------------------------------------------------------
54      INTEGER, INTENT(in) ::   kt    ! number of iteration
55      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
56      REAL(wp) ::   zsal
57      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
58      !!-------------------------------------------------------------------
59      IF( nn_timing == 1 )  CALL timing_start('limupdate2')
60
61      IF( kt == nit000 .AND. lwp ) THEN
62         WRITE(numout,*)''
63         WRITE(numout,*)' lim_update2 '
64         WRITE(numout,*)' ~~~~~~~~~~~ '
65      ENDIF
66
67      ! conservation test
68      IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
69
70      !----------------------------------------------------------------------
71      ! Constrain the thickness of the smallest category above himin
72      !----------------------------------------------------------------------
73!$OMP PARALLEL
74!$OMP DO schedule(static) private(jj,ji,rswitch)
75      DO jj = 1, jpj 
76         DO ji = 1, jpi
77            rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes
78            ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch
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               oa_i(ji,jj,1) = oa_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin
82            ENDIF
83         END DO
84      END DO
85     
86      !-----------------------------------------------------
87      ! ice concentration should not exceed amax
88      !-----------------------------------------------------
89!$OMP DO schedule(static) private(jj, ji)
90      DO jj = 1, jpj
91         DO ji = 1, jpi
92            at_i(ji,jj) = 0._wp
93         END DO
94      END DO
95      DO jl = 1, jpl
96!$OMP DO schedule(static) private(jj, ji)
97         DO jj = 1, jpj
98            DO ji = 1, jpi
99               at_i(ji,jj) = a_i(ji,jj,jl) + at_i(ji,jj)
100            END DO
101         END DO
102      END DO
103
104      DO jl  = 1, jpl
105!$OMP DO schedule(static) private(jj, ji)
106         DO jj = 1, jpj
107            DO ji = 1, jpi
108               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
109                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
110                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
111               ENDIF
112            END DO
113         END DO
114      END DO
115!$OMP END PARALLEL
116
117      !---------------------
118      ! Ice salinity
119      !---------------------
120      IF (  nn_icesal == 2  ) THEN
121         DO jl = 1, jpl
122!$OMP PARALLEL DO schedule(static) private(jj,ji,zsal,rswitch)
123            DO jj = 1, jpj 
124               DO ji = 1, jpi
125                  zsal            = smv_i(ji,jj,jl)
126                  ! salinity stays in bounds
127                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
128                  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) )
129                  ! associated salt flux
130                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
131               END DO
132            END DO
133         END DO
134      ENDIF
135
136      !----------------------------------------------------
137      ! Rebin categories with thickness out of bounds
138      !----------------------------------------------------
139      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )
140
141      !-----------------
142      ! zap small values
143      !-----------------
144      CALL lim_var_zapsmall
145
146      !------------------------------------------------------------------------------
147      ! Corrections to avoid wrong values                                        |
148      !------------------------------------------------------------------------------
149      ! Ice drift
150      !------------
151!$OMP PARALLEL DO schedule(static) private(jj, ji)
152      DO jj = 2, jpjm1
153         DO ji = 2, jpim1
154            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
155               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
156               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
157               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
158               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
159            ENDIF
160         END DO
161      END DO
162      !lateral boundary conditions
163      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
164      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
165      !mask velocities
166!$OMP PARALLEL
167!$OMP DO schedule(static) private(jj, ji)
168      DO jj = 1, jpj
169         DO ji = 1, jpi
170            u_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1)
171            v_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
172         END DO
173      END DO
174 
175      ! -------------------------------------------------
176      ! Diagnostics
177      ! -------------------------------------------------
178      DO jl  = 1, jpl
179!$OMP DO schedule(static) private(jj, ji)
180         DO jj = 1, jpj
181            DO ji = 1, jpi
182               oa_i(ji,jj,jl) = oa_i(ji,jj,jl) + a_i(ji,jj,jl) * rdt_ice / rday   ! ice natural aging
183               afx_thd(ji,jj) = afx_thd(ji,jj) + ( a_i(ji,jj,jl) - a_i_b(ji,jj,jl) ) * r1_rdtice
184            END DO
185         END DO
186      END DO
187      afx_tot = afx_thd + afx_dyn
188
189!$OMP DO schedule(static) private(jj, ji)
190      DO jj = 1, jpj
191         DO ji = 1, jpi           
192            ! heat content variation (W.m-2)
193            diag_heat(ji,jj) = diag_heat(ji,jj) -  &
194               &               ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
195               &                 SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
196               &               ) * r1_rdtice   
197            ! salt, volume
198            diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice
199            diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice
200            diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice
201         END DO
202      END DO
203!$OMP END PARALLEL
204
205      ! conservation test
206      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
207
208      ! control prints
209      IF( ln_limctl )    CALL lim_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
210      IF( ln_ctl )       CALL lim_prt3D( 'limupdate2' )
211   
212      IF( nn_timing == 1 )  CALL timing_stop('limupdate2')
213
214   END SUBROUTINE lim_update2
215
216#else
217   !!----------------------------------------------------------------------
218   !!   Default option         Empty Module               No sea-ice model
219   !!----------------------------------------------------------------------
220CONTAINS
221   SUBROUTINE lim_update2     ! Empty routine
222   END SUBROUTINE lim_update2
223
224#endif
225
226END MODULE limupdate2
Note: See TracBrowser for help on using the repository browser.