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/UKMO/dev_10448_WAD_SBC_BUGFIX/src/OCE/DOM – NEMO

source: NEMO/branches/UKMO/dev_10448_WAD_SBC_BUGFIX/src/OCE/DOM/istate.F90 @ 10456

Last change on this file since 10456 was 10456, checked in by deazer, 5 years ago

Added option to taper sbc fluxes near very shallow water when using WAD
Corrected some IO bugs in dia25h, diatmb for WAD case.
User has control of the tapering. At what depth to start it, and at what fraction to start
the tanh tapering. At the WAD limit SBC is turned off completely.
Dry cells do not have any communication with the atmosphere
To DO: Documentation update.
Although not all sette tests are passed (AGRIF etc.)
it does no worse than the trunk at the revision the branch is made
  • Property svn:keywords set to Id
File size: 8.5 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_tsd_init ) THEN               
99            CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
100            !
101            sshb(:,:)   = 0._wp               ! set the ocean at rest
102            IF(ll_wd) then
103               sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD
104               !
105               ! Apply minimum wetdepth criterion
106               !
107               DO jj = 1,jpj
108                  DO ji = 1,jpi
109                     IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN
110                        sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
111                     ENDIF
112                  END DO
113               END DO
114            ENDIF
115            ub  (:,:,:) = 0._wp
116            vb  (:,:,:) = 0._wp 
117            !
118         ELSE                                 ! user defined initial T and S
119            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
120         ENDIF
121         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
122         sshn (:,:)     = sshb(:,:)   
123         un   (:,:,:)   = ub  (:,:,:)
124         vn   (:,:,:)   = vb  (:,:,:)
125         hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
126         CALL div_hor( 0 )                    ! compute interior hdivn value 
127!!gm                                    hdivn(:,:,:) = 0._wp
128
129!!gm POTENTIAL BUG :
130!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
131!!             as well as gdept and gdepw....   !!!!!
132!!      ===>>>>   probably a call to domvvl initialisation here....
133
134
135         !
136!!gm to be moved in usrdef of C1D case
137!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
138!            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
139!            CALL dta_uvd( nit000, zuvd )
140!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
141!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
142!            DEALLOCATE( zuvd )
143!         ENDIF
144         !
145!!gm This is to be changed !!!!
146!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
147!         IF( .NOT.ln_linssh ) THEN
148!            DO jk = 1, jpk
149!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
150!            END DO
151!         ENDIF
152!!gm
153         !
154      ENDIF 
155      !
156      ! Initialize "now" and "before" barotropic velocities:
157      ! Do it whatever the free surface method, these arrays being eventually used
158      !
159      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
160      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
161      !
162!!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
163      DO jk = 1, jpkm1
164         DO jj = 1, jpj
165            DO ji = 1, jpi
166               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
167               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
168               !
169               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
170               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
171            END DO
172         END DO
173      END DO
174      !
175      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
176      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
177      !
178      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
179      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
180      !
181   END SUBROUTINE istate_init
182
183   !!======================================================================
184END MODULE istate
Note: See TracBrowser for help on using the repository browser.