source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 7698

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

update trunk with OpenMP parallelization

  • Property svn:keywords set to Id
File size: 10.0 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 iscplrst        ! ice sheet coupling
31   USE wet_dry         ! wetting and drying (needed for wad_istate)
32   USE usrdef_istate   ! User defined initial state
33   !
34   USE in_out_manager  ! I/O manager
35   USE iom             ! I/O library
36   USE lib_mpp         ! MPP library
37   USE restart         ! restart
38   USE wrk_nemo        ! Memory allocation
39   USE timing          ! Timing
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   istate_init   ! routine called by step.F90
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE istate_init
56      !!----------------------------------------------------------------------
57      !!                   ***  ROUTINE istate_init  ***
58      !!
59      !! ** Purpose :   Initialization of the dynamics and tracer fields.
60      !!----------------------------------------------------------------------
61      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
62      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
63      !!----------------------------------------------------------------------
64      !
65      IF( nn_timing == 1 )   CALL timing_start('istate_init')
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!$OMP PARALLEL
78      DO jn = 1, jpts
79!$OMP DO schedule(static) private(jk, jj, ji)
80         DO jk = 1, jpk
81            DO jj = 1, jpj
82               DO ji = 1, jpi
83                  tsa  (ji,jj,jk,jn) = 0._wp                                       ! set one for all to 0 at level jpk
84                  rab_b(ji,jj,jk,jn) = 0._wp   ;   rab_n(ji,jj,jk,jn) = 0._wp      ! set one for all to 0 at level jpk
85               END DO
86            END DO
87         END DO
88      END DO
89!$OMP DO schedule(static) private(jk, jj, ji)
90      DO jk = 1, jpk
91         DO jj = 1, jpj
92            DO ji = 1, jpi
93               rhd  (ji,jj,jk  ) = 0._wp   ;   rhop (ji,jj,jk  ) = 0._wp      ! set one for all to 0 at level jpk
94               rn2b (ji,jj,jk  ) = 0._wp   ;   rn2  (ji,jj,jk  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
95            END DO
96         END DO
97      END DO
98!$OMP END PARALLEL
99
100      IF( ln_rstart ) THEN                    ! Restart from a file
101         !                                    ! -------------------
102         CALL rst_read                        ! Read the restart file
103         IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry
104         CALL day_init                        ! model calendar (using both namelist and restart infos)
105         !
106      ELSE                                    ! Start from rest
107         !                                    ! ---------------
108         numror = 0                           ! define numror = 0 -> no restart file to read
109         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
110         CALL day_init                        ! model calendar (using both namelist and restart infos)
111         !                                    ! Initialization of ocean to zero
112         !
113         IF( ln_tsd_init ) THEN               
114            CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
115            !
116!$OMP PARALLEL
117!$OMP DO schedule(static) private(jj, ji)
118            DO jj = 1, jpj
119               DO ji = 1, jpi
120                  sshb (ji,jj)   = 0._wp      ! set the ocean at rest
121               END DO
122            END DO
123!$OMP END DO NOWAIT
124!$OMP DO schedule(static) private(jk, jj, ji)
125            DO jk = 1, jpk
126               DO jj = 1, jpj
127                  DO ji = 1, jpi
128                     ub  (ji,jj,jk) = 0._wp
129                     vb  (ji,jj,jk) = 0._wp 
130                  END DO
131               END DO
132            END DO
133!$OMP END PARALLEL
134            !
135         ELSE                                 ! user defined initial T and S
136            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
137         ENDIF
138!$OMP PARALLEL
139         DO jn = 1, jpts
140!$OMP DO schedule(static) private(jk, jj, ji)
141            DO jk = 1, jpk
142               DO jj = 1, jpj
143                  DO ji = 1, jpi
144                     tsn  (ji,jj,jk,jn) = tsb (ji,jj,jk,jn)       ! set now values from to before ones
145                  END DO
146               END DO
147            END DO
148         END DO
149!$OMP DO schedule(static) private(jk, jj, ji)
150         DO jk = 1, jpk
151            DO jj = 1, jpj
152               DO ji = 1, jpi
153                  un   (ji,jj,jk)   = ub  (ji,jj,jk)
154                  vn   (ji,jj,jk)   = vb  (ji,jj,jk)
155               END DO
156            END DO
157         END DO
158!$OMP DO schedule(static) private(jj, ji)
159         DO jj = 1, jpj
160            DO ji = 1, jpi
161               sshn (ji,jj)     = sshb(ji,jj)   
162               hdivn(ji,jj,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
163            END DO
164         END DO
165!$OMP END PARALLEL
166         CALL div_hor( 0 )                    ! compute interior hdivn value 
167!!gm                                    hdivn(:,:,:) = 0._wp
168
169!!gm POTENTIAL BUG :
170!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
171!!             as well as gdept and gdepw....   !!!!!
172!!      ===>>>>   probably a call to domvvl initialisation here....
173
174
175         !
176!!gm to be moved in usrdef of C1D case
177!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
178!            CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
179!            CALL dta_uvd( nit000, zuvd )
180!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
181!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
182!            CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
183!         ENDIF
184         !
185!!gm This is to be changed !!!!
186!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
187!         IF( .NOT.ln_linssh ) THEN
188!            DO jk = 1, jpk
189!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
190!            END DO
191!         ENDIF
192!!gm
193         !
194      ENDIF 
195      !
196      ! Initialize "now" and "before" barotropic velocities:
197      ! Do it whatever the free surface method, these arrays being eventually used
198      !
199!$OMP PARALLEL
200!$OMP DO schedule(static) private(jj, ji)
201      DO jj = 1, jpj
202         DO ji = 1, jpi
203            un_b(ji,jj) = 0._wp   ;   vn_b(ji,jj) = 0._wp
204            ub_b(ji,jj) = 0._wp   ;   vb_b(ji,jj) = 0._wp
205         END DO
206      END DO
207      !
208!!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
209      DO jk = 1, jpkm1
210!$OMP DO schedule(static) private(jj, ji)
211         DO jj = 1, jpj
212            DO ji = 1, jpi
213               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
214               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
215               !
216               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
217               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
218            END DO
219         END DO
220      END DO
221      !
222!$OMP DO schedule(static) private(jj, ji)
223      DO jj = 1, jpj
224         DO ji = 1, jpi
225            un_b(ji,jj) = un_b(ji,jj) * r1_hu_n(ji,jj)
226            vn_b(ji,jj) = vn_b(ji,jj) * r1_hv_n(ji,jj)
227            !
228            ub_b(ji,jj) = ub_b(ji,jj) * r1_hu_b(ji,jj)
229            vb_b(ji,jj) = vb_b(ji,jj) * r1_hv_b(ji,jj)
230         END DO
231      END DO
232!$OMP END PARALLEL
233      !
234      IF( nn_timing == 1 )   CALL timing_stop('istate_init')
235      !
236   END SUBROUTINE istate_init
237
238   !!======================================================================
239END MODULE istate
Note: See TracBrowser for help on using the repository browser.