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/releases/r4.0/r4.0-HEAD/src/OCE/DOM – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/istate.F90 @ 15458

Last change on this file since 15458 was 15458, checked in by clem, 3 years ago

4.0-HEAD: recalculate vertical scale factors in istate.F90 for wetting and drying and user defined configurations (related to change of ssh). Should solve ticket #2736 for 4.0. Not sure we also need to change the trunk?
  • Property svn:keywords set to Id
File size: 9.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 dtatsd         ! data temperature and salinity   (dta_tsd routine)
27   USE dtauvd         ! data: U & V current             (dta_uvd routine)
28   USE domvvl          ! varying vertical mesh
29   USE iscplrst        ! ice sheet coupling
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   !!----------------------------------------------------------------------
46   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
47   !! $Id$
48   !! Software governed by the CeCILL license (see ./LICENSE)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE istate_init
53      !!----------------------------------------------------------------------
54      !!                   ***  ROUTINE istate_init  ***
55      !!
56      !! ** Purpose :   Initialization of the dynamics and tracer fields.
57      !!----------------------------------------------------------------------
58      INTEGER ::   ji, jj, jk   ! dummy loop indices
59!!gm see comment further down
60      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
61!!gm end
62      REAL(wp), DIMENSION(jpi,jpj) ::   z2d    ! temporary array
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               !
115               IF( .NOT.ln_linssh ) THEN      ! adjust initial vertical scale factors
116                  sshn(:,:) = sshb(:,:)
117                  z2d (:,:) = 1._wp + sshn(:,:) * ssmask(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) )
118                  DO jk = 1, jpkm1               
119                     DO jj = 1, jpj
120                        DO ji = 1, jpi
121                           IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN
122                              e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( z2d(ji,jj) * tmask(ji,jj,jk) - ( tmask(ji,jj,jk) - 1._wp ) )
123                              e3t_b(ji,jj,jk) = e3t_n(ji,jj,jk)
124                              e3t_a(ji,jj,jk) = e3t_n(ji,jj,jk)
125                           ENDIF
126                        END DO
127                     END DO
128                  END DO
129                  CALL dom_vvl_zgr
130               ENDIF
131               !
132            ENDIF
133            ub  (:,:,:) = 0._wp
134            vb  (:,:,:) = 0._wp
135            !
136         ELSE                                 ! user defined initial T and S
137            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
138            !
139            IF( .NOT.ln_linssh ) THEN         ! adjust initial vertical scale factors
140               sshn(:,:) = sshb(:,:)
141               z2d (:,:) = 1._wp + sshn(:,:) * ssmask(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) )
142               DO jk = 1, jpkm1               
143                  DO jj = 1, jpj
144                     DO ji = 1, jpi
145                        e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( z2d(ji,jj) * tmask(ji,jj,jk) - ( tmask(ji,jj,jk) - 1._wp ) )
146                        e3t_b(ji,jj,jk) = e3t_n(ji,jj,jk)
147                        e3t_a(ji,jj,jk) = e3t_n(ji,jj,jk)
148                     END DO
149                  END DO
150               END DO
151               CALL dom_vvl_zgr
152            ENDIF
153            !
154         ENDIF
155         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
156         sshn (:,:)     = sshb(:,:)   
157         un   (:,:,:)   = ub  (:,:,:)
158         vn   (:,:,:)   = vb  (:,:,:)
159
160!!gm POTENTIAL BUG :
161!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
162!!             as well as gdept and gdepw....   !!!!!
163!!      ===>>>>   probably a call to domvvl initialisation here....
164
165
166         !
167!!gm to be moved in usrdef of C1D case
168!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
169!            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
170!            CALL dta_uvd( nit000, zuvd )
171!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
172!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
173!            DEALLOCATE( zuvd )
174!         ENDIF
175         !
176!!gm This is to be changed !!!!
177!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
178!         IF( .NOT.ln_linssh ) THEN
179!            DO jk = 1, jpk
180!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
181!            END DO
182!         ENDIF
183!!gm
184         !
185      ENDIF 
186      !
187      ! Initialize "now" and "before" barotropic velocities:
188      ! Do it whatever the free surface method, these arrays being eventually used
189      !
190      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
191      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
192      !
193!!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
194      DO jk = 1, jpkm1
195         DO jj = 1, jpj
196            DO ji = 1, jpi
197               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
198               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
199               !
200               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
201               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
202            END DO
203         END DO
204      END DO
205      !
206      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
207      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
208      !
209      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
210      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
211      !
212   END SUBROUTINE istate_init
213
214   !!======================================================================
215END MODULE istate
Note: See TracBrowser for help on using the repository browser.