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.
istate.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/istate.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: 8.7 KB
Line 
1MODULE istate
2   !!======================================================================
3   !!                     ***  MODULE  istate  ***
4   !! Ocean state   :  initial state setting
5   !!=====================================================================
6   !! History :  OPA  !  1989-12  (P. Andrich)  Original code
7   !!            5.0  !  1991-11  (G. Madec)  rewritting
8   !!            6.0  !  1996-01  (G. Madec)  terrain following coordinates
9   !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_eel
10   !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_uvg
11   !!   NEMO     1.0  !  2003-08  (G. Madec, C. Talandier)  F90: Free form, modules + EEL R5
12   !!             -   !  2004-05  (A. Koch-Larrouy)  istate_gyre
13   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
14   !!            3.3  !  2010-10  (C. Ethe) merge TRC-TRA
15   !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn
16   !!            3.7  !  2016-04  (S. Flavoni) introduce user defined initial state
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   istate_init   : initial state setting
21   !!   istate_uvg    : initial velocity in geostropic balance
22   !!----------------------------------------------------------------------
23   USE oce            ! ocean dynamics and active tracers
24   USE dom_oce        ! ocean space and time domain
25   USE daymod         ! calendar
26   USE divhor         ! horizontal divergence            (div_hor routine)
27   USE dtatsd         ! data temperature and salinity   (dta_tsd routine)
28   USE dtauvd         ! data: U & V current             (dta_uvd routine)
29   USE domvvl          ! varying vertical mesh
30   USE wet_dry         ! wetting and drying (needed for wad_istate)
31   USE usrdef_istate   ! User defined initial state
32   !
33   USE in_out_manager  ! I/O manager
34   USE iom             ! I/O library
35   USE lib_mpp         ! MPP library
36   USE restart         ! restart
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   istate_init   ! routine called by step.F90
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45#  include "do_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
48   !! $Id$
49   !! Software governed by the CeCILL license (see ./LICENSE)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE istate_init( Kbb, Kmm, Kaa )
54      !!----------------------------------------------------------------------
55      !!                   ***  ROUTINE istate_init  ***
56      !!
57      !! ** Purpose :   Initialization of the dynamics and tracer fields.
58      !!----------------------------------------------------------------------
59      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices
60      !
61      INTEGER ::   ji, jj, jk   ! dummy loop indices
62!!gm see comment further down
63      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
64!!gm end
65      !!----------------------------------------------------------------------
66      !
67      IF(lwp) WRITE(numout,*)
68      IF(lwp) WRITE(numout,*) 'istate_init : Initialization of the dynamics and tracers'
69      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
70
71!!gm  Why not include in the first call of dta_tsd ? 
72!!gm  probably associated with the use of internal damping...
73                     CALL dta_tsd_init        ! Initialisation of T & S input data
74!!gm to be moved in usrdef of C1D case
75!      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
76!!gm
77
78      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
79      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
80      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
81      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
82#if defined key_agrif
83      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization
84      vv   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization   
85#endif
86
87      IF( ln_rstart ) THEN                    ! Restart from a file
88         !                                    ! -------------------
89         CALL rst_read( Kbb, Kmm )            ! Read the restart file
90         CALL day_init                        ! model calendar (using both namelist and restart infos)
91         !
92      ELSE                                    ! Start from rest
93         !                                    ! ---------------
94         numror = 0                           ! define numror = 0 -> no restart file to read
95         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
96         CALL day_init                        ! model calendar (using both namelist and restart infos)
97         !                                    ! Initialization of ocean to zero
98         !
99         IF( ln_tsd_init ) THEN               
100            CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000
101            !
102            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest
103            IF( ll_wd ) THEN
104               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD
105               !
106               ! Apply minimum wetdepth criterion
107               !
108               DO_2D_11_11
109                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
110                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
111                  ENDIF
112               END_2D
113            ENDIF
114            uu  (:,:,:,Kbb) = 0._wp
115            vv  (:,:,:,Kbb) = 0._wp 
116            !
117         ELSE                                 ! user defined initial T and S
118            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
119         ENDIF
120         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
121         ssh (:,:,Kmm)     = ssh(:,:,Kbb)   
122         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
123         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
124         hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
125         CALL div_hor( 0, Kbb, Kmm )         ! compute interior hdiv value 
126!!gm                                    hdiv(:,:,:) = 0._wp
127
128!!gm POTENTIAL BUG :
129!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
130!!             as well as gdept and gdepw....   !!!!!
131!!      ===>>>>   probably a call to domvvl initialisation here....
132
133
134         !
135!!gm to be moved in usrdef of C1D case
136!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
137!            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
138!            CALL dta_uvd( nit000, zuvd )
139!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
140!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
141!            DEALLOCATE( zuvd )
142!         ENDIF
143         !
144!!gm This is to be changed !!!!
145!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here
146!         IF( .NOT.ln_linssh ) THEN
147!            DO jk = 1, jpk
148!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
149!            END DO
150!         ENDIF
151!!gm
152         !
153      ENDIF 
154      !
155      ! Initialize "now" and "before" barotropic velocities:
156      ! Do it whatever the free surface method, these arrays being eventually used
157      !
158      uu_b(:,:,Kmm) = 0._wp   ;   vv_b(:,:,Kmm) = 0._wp
159      uu_b(:,:,Kbb) = 0._wp   ;   vv_b(:,:,Kbb) = 0._wp
160      !
161!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
162      DO_3D_11_11( 1, jpkm1 )
163         uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
164         vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
165         !
166         uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
167         vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
168      END_3D
169      !
170      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
171      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
172      !
173      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb)
174      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb)
175      !
176   END SUBROUTINE istate_init
177
178   !!======================================================================
179END MODULE istate
Note: See TracBrowser for help on using the repository browser.