source: NEMO/branches/UKMO/NEMO_4.0_surge/src/OCE/DOM/istate.F90 @ 11180

Last change on this file since 11180 was 11180, checked in by clne, 22 months ago

Initial commit of code for 2d (surge) work in NEMO4.
This is aiming to replicate the 3.6 version in branches/UKMO/dev_r5518_Surge_Modelling

File size: 8.9 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
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   istate_init   ! routine called by step.F90
43
44   !! * Substitutions
45#  include "vectopt_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
54      !!----------------------------------------------------------------------
55      !!                   ***  ROUTINE istate_init  ***
56      !!
57      !! ** Purpose :   Initialization of the dynamics and tracer fields.
58      !!----------------------------------------------------------------------
59      INTEGER ::   ji, jj, jk   ! dummy loop indices
60!!gm see comment further down
61      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
62!!gm end
63      !!----------------------------------------------------------------------
64      !
65      IF(lwp) WRITE(numout,*)
66      IF(lwp) WRITE(numout,*) 'istate_init : Initialization of the dynamics and tracers'
67      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
68
69!!gm  Why not include in the first call of dta_tsd ? 
70!!gm  probably associated with the use of internal damping...
71                     CALL dta_tsd_init        ! Initialisation of T & S input data
72!!gm to be moved in usrdef of C1D case
73!      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
74!!gm
75
76      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
77      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
78      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
79      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
80#if defined key_agrif
81      ua   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization
82      va   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization   
83#endif
84
85      IF( ln_rstart ) THEN                    ! Restart from a file
86         !                                    ! -------------------
87         CALL rst_read                        ! Read the restart file
88         IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry
89         CALL day_init                        ! model calendar (using both namelist and restart infos)
90         !
91      ELSE                                    ! Start from rest
92         !                                    ! ---------------
93         numror = 0                           ! define numror = 0 -> no restart file to read
94         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
95         CALL day_init                        ! model calendar (using both namelist and restart infos)
96         !                                    ! Initialization of ocean to zero
97         !
98         IF (ln_2d) THEN
99            IF(lwp) WRITE(numout,*) 'istate_init : 2D case, setting tracers to contants and ocean at rest' 
100            tsb(:,:,:,:)= 15._wp               ! No tracers, but can't set salinity to 0 otherwise it triggers crash message
101            sshb(:,:)   = 0._wp                ! set the ocean at rest
102            ub  (:,:,:) = 0._wp               
103            vb  (:,:,:) = 0._wp           
104         ELSE IF( ln_tsd_init ) THEN   
105            CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
106            !
107            sshb(:,:)   = 0._wp               ! set the ocean at rest
108            IF( ll_wd ) THEN
109               sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD
110               !
111               ! Apply minimum wetdepth criterion
112               !
113               DO jj = 1,jpj
114                  DO ji = 1,jpi
115                     IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN
116                        sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
117                     ENDIF
118                  END DO
119               END DO
120            ENDIF
121            ub  (:,:,:) = 0._wp
122            vb  (:,:,:) = 0._wp 
123            !
124         ELSE                                 ! user defined initial T and S
125            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
126         ENDIF
127         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
128         sshn (:,:)     = sshb(:,:)   
129         un   (:,:,:)   = ub  (:,:,:)
130         vn   (:,:,:)   = vb  (:,:,:)
131         hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
132         CALL div_hor( 0 )                    ! compute interior hdivn value 
133!!gm                                    hdivn(:,:,:) = 0._wp
134
135!!gm POTENTIAL BUG :
136!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
137!!             as well as gdept and gdepw....   !!!!!
138!!      ===>>>>   probably a call to domvvl initialisation here....
139
140
141         !
142!!gm to be moved in usrdef of C1D case
143!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
144!            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
145!            CALL dta_uvd( nit000, zuvd )
146!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
147!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
148!            DEALLOCATE( zuvd )
149!         ENDIF
150         !
151!!gm This is to be changed !!!!
152!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
153!         IF( .NOT.ln_linssh ) THEN
154!            DO jk = 1, jpk
155!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
156!            END DO
157!         ENDIF
158!!gm
159         !
160      ENDIF 
161      !
162      ! Initialize "now" and "before" barotropic velocities:
163      ! Do it whatever the free surface method, these arrays being eventually used
164      !
165      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
166      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
167      !
168!!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
169      DO jk = 1, jpkm1
170         DO jj = 1, jpj
171            DO ji = 1, jpi
172               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
173               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
174               !
175               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
176               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
177            END DO
178         END DO
179      END DO
180      !
181      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
182      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
183      !
184      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
185      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
186      !
187   END SUBROUTINE istate_init
188
189   !!======================================================================
190END MODULE istate
Note: See TracBrowser for help on using the repository browser.