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.
icedyn.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 14.3 KB
Line 
1MODULE icedyn
2   !!======================================================================
3   !!                     ***  MODULE  icedyn  ***
4   !!   Sea-Ice dynamics : master routine for sea ice dynamics
5   !!======================================================================
6   !! history :  4.0  ! 2018  (C. Rousset)  original code SI3 [aka Sea Ice cube]
7   !!----------------------------------------------------------------------
8#if defined key_si3
9   !!----------------------------------------------------------------------
10   !!   'key_si3'                                       SI3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!   ice_dyn       : dynamics of sea ice
13   !!   ice_dyn_init  : initialization and namelist read
14   !!----------------------------------------------------------------------
15   USE phycst         ! physical constants
16   USE dom_oce        ! ocean space and time domain
17   USE ice            ! sea-ice: variables
18   USE icedyn_rhg     ! sea-ice: rheology
19   USE icedyn_adv     ! sea-ice: advection
20   USE icedyn_rdgrft  ! sea-ice: ridging/rafting
21   USE icecor         ! sea-ice: corrections
22   USE icevar         ! sea-ice: operations
23   USE icectl         ! sea-ice: control prints
24   !
25   USE in_out_manager ! I/O manager
26   USE iom            ! I/O manager library
27   USE lib_mpp        ! MPP library
28   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
29   USE lbclnk         ! lateral boundary conditions (or mpp links)
30   USE timing         ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   ice_dyn        ! called by icestp.F90
36   PUBLIC   ice_dyn_init   ! called by icestp.F90
37   
38   INTEGER ::              nice_dyn   ! choice of the type of dynamics
39   !                                        ! associated indices:
40   INTEGER, PARAMETER ::   np_dynALL     = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
41   INTEGER, PARAMETER ::   np_dynRHGADV  = 2   ! pure dynamics                   (rheology + advection)
42   INTEGER, PARAMETER ::   np_dynADV1D   = 3   ! only advection 1D - test case from Schar & Smolarkiewicz 1996
43   INTEGER, PARAMETER ::   np_dynADV2D   = 4   ! only advection 2D w prescribed vel.(rn_uvice + advection)
44   !
45   ! ** namelist (namdyn) **
46   LOGICAL  ::   ln_dynALL        ! full ice dynamics                      (rheology + advection + ridging/rafting + correction)
47   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections         (rheology + advection)
48   LOGICAL  ::   ln_dynADV1D      ! only advection in 1D w ice convergence (test case from Schar & Smolarkiewicz 1996)
49   LOGICAL  ::   ln_dynADV2D      ! only advection in 2D w prescribed vel. (rn_uvice + advection)
50   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV1D & np_dynADV2D)
51   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV1D & np_dynADV2D)
52   
53   !! * Substitutions
54#  include "vectopt_loop_substitute.h90"
55#  include "do_loop_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (./LICENSE)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE ice_dyn( kt, Kmm )
64      !!-------------------------------------------------------------------
65      !!               ***  ROUTINE ice_dyn  ***
66      !!               
67      !! ** Purpose :   this routine manages sea ice dynamics
68      !!
69      !! ** Action : - calculation of friction in case of landfast ice
70      !!             - call ice_dyn_rhg    = rheology
71      !!             - call ice_dyn_adv    = advection
72      !!             - call ice_dyn_rdgrft = ridging/rafting
73      !!             - call ice_cor        = corrections if fields are out of bounds
74      !!--------------------------------------------------------------------
75      INTEGER, INTENT(in) ::   kt     ! ice time step
76      INTEGER, INTENT(in) ::   Kmm    ! ocean time level index
77      !!
78      INTEGER  ::   ji, jj        ! dummy loop indices
79      REAL(wp) ::   zcoefu, zcoefv
80      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i
81      !!--------------------------------------------------------------------
82      !
83      ! controls
84      IF( ln_timing )   CALL timing_start('icedyn')
85      !
86      IF( kt == nit000 .AND. lwp ) THEN
87         WRITE(numout,*)
88         WRITE(numout,*)'ice_dyn: sea-ice dynamics'
89         WRITE(numout,*)'~~~~~~~'
90      ENDIF
91      !                     
92      ! retrieve thickness from volume for landfast param. and UMx advection scheme
93      WHERE( a_i(:,:,:) >= epsi20 )
94         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
95         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
96      ELSEWHERE
97         h_i(:,:,:) = 0._wp
98         h_s(:,:,:) = 0._wp
99      END WHERE
100      !
101      WHERE( a_ip(:,:,:) >= epsi20 )
102         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
103      ELSEWHERE
104         h_ip(:,:,:) = 0._wp
105      END WHERE
106      !
107      !
108      SELECT CASE( nice_dyn )          !-- Set which dynamics is running
109
110      CASE ( np_dynALL )           !==  all dynamical processes  ==!
111         !
112         CALL ice_dyn_rhg   ( kt, Kmm )                                     ! -- rheology 
113         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
114         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting
115         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections
116         !
117      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
118         !
119         CALL ice_dyn_rhg   ( kt, Kmm )                                     ! -- rheology 
120         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
121         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting)
122         CALL ice_var_zapsmall                                              ! -- zap small areas
123         !
124      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
125         !
126         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
127         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
128         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s
129         DO_2D_11_11
130            zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
131            zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
132            u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
133            v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
134         END_2D
135         ! ---
136         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
137         !
138      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities)
139         !
140         u_ice(:,:) = rn_uice * umask(:,:,1)
141         v_ice(:,:) = rn_vice * vmask(:,:,1)
142         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
143         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
144         ! ---
145         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
146
147      END SELECT
148      !
149      !
150      ! diagnostics: divergence at T points
151      IF( iom_use('icediv') ) THEN
152         !
153         SELECT CASE( nice_dyn )
154
155         CASE ( np_dynADV1D , np_dynADV2D )
156
157            ALLOCATE( zdivu_i(jpi,jpj) )
158            DO_2D_00_00
159               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
160                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
161            END_2D
162            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
163            ! output
164            CALL iom_put( 'icediv' , zdivu_i )
165
166            DEALLOCATE( zdivu_i )
167
168         END SELECT
169         !
170      ENDIF
171      !
172      ! controls
173      IF( ln_timing )   CALL timing_stop ('icedyn')
174      !
175   END SUBROUTINE ice_dyn
176
177
178   SUBROUTINE Hpiling
179      !!-------------------------------------------------------------------
180      !!                  ***  ROUTINE Hpiling  ***
181      !!
182      !! ** Purpose : Simple conservative piling comparable with 1-cat models
183      !!
184      !! ** Method  : pile-up ice when no ridging/rafting
185      !!
186      !! ** input   : a_i
187      !!-------------------------------------------------------------------
188      INTEGER ::   jl         ! dummy loop indices
189      !!-------------------------------------------------------------------
190      ! controls
191      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
192      !
193      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
194      DO jl = 1, jpl
195         WHERE( at_i(:,:) > epsi20 )
196            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
197         END WHERE
198      END DO
199      ! controls
200      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
201      !
202   END SUBROUTINE Hpiling
203
204
205   SUBROUTINE ice_dyn_init
206      !!-------------------------------------------------------------------
207      !!                  ***  ROUTINE ice_dyn_init  ***
208      !!
209      !! ** Purpose : Physical constants and parameters linked to the ice
210      !!      dynamics
211      !!
212      !! ** Method  :  Read the namdyn namelist and check the ice-dynamic
213      !!       parameter values called at the first timestep (nit000)
214      !!
215      !! ** input   :   Namelist namdyn
216      !!-------------------------------------------------------------------
217      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
218      !!
219      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  &
220         &             rn_ishlat ,                                                           &
221         &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
222      !!-------------------------------------------------------------------
223      !
224      READ  ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
225901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn in reference namelist' )
226      READ  ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
227902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn in configuration namelist' )
228      IF(lwm) WRITE( numoni, namdyn )
229      !
230      IF(lwp) THEN                     ! control print
231         WRITE(numout,*)
232         WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
233         WRITE(numout,*) '~~~~~~~~~~~~'
234         WRITE(numout,*) '   Namelist namdyn:'
235         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL
236         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV
237         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D
238         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D
239         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')'
240         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
241         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
242         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
243         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
244         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax
245         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile
246         WRITE(numout,*)
247      ENDIF
248      !                             !== set the choice of ice dynamics ==!
249      ioptio = 0 
250      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
251      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF
252      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
253      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF
254      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996
255      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF
256      !      !--- advection 2D only with prescribed ice velocities (from namelist)
257      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF
258      !
259      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
260      !
261      !                                      !--- Lateral boundary conditions
262      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip'
263      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip'
264      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip'
265      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
266      ENDIF
267      !                                      !--- Landfast ice
268      IF( .NOT.ln_landfast_L16 )   tau_icebfr(:,:) = 0._wp
269      !
270      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
271      CALL ice_dyn_rhg_init             ! set ice rheology parameters
272      CALL ice_dyn_adv_init             ! set ice advection parameters
273      !
274   END SUBROUTINE ice_dyn_init
275
276#else
277   !!----------------------------------------------------------------------
278   !!   Default option         Empty module           NO SI3 sea-ice model
279   !!----------------------------------------------------------------------
280#endif 
281
282   !!======================================================================
283END MODULE icedyn
Note: See TracBrowser for help on using the repository browser.