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 branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/DOM – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/DOM/istate.F90 @ 8568

Last change on this file since 8568 was 8568, checked in by gm, 7 years ago

#1911 (ENHANCE-09): PART I.2 - _NONE option + remove zts + see associated wiki page

  • Property svn:keywords set to Id
File size: 8.0 KB
RevLine 
[3]1MODULE istate
2   !!======================================================================
3   !!                     ***  MODULE  istate  ***
4   !! Ocean state   :  initial state setting
5   !!=====================================================================
[2104]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
[3294]15   !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn
[7646]16   !!            3.7  !  2016-04  (S. Flavoni) introduce user defined initial state
[508]17   !!----------------------------------------------------------------------
[3]18
19   !!----------------------------------------------------------------------
20   !!   istate_init   : initial state setting
21   !!   istate_uvg    : initial velocity in geostropic balance
22   !!----------------------------------------------------------------------
[7646]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)
[593]29   USE domvvl          ! varying vertical mesh
[6140]30   USE iscplrst        ! ice sheet coupling
[7646]31   USE wet_dry         ! wetting and drying (needed for wad_istate)
32   USE usrdef_istate   ! User defined initial state
[4990]33   !
34   USE in_out_manager  ! I/O manager
35   USE iom             ! I/O library
[2715]36   USE lib_mpp         ! MPP library
[3680]37   USE restart         ! restart
[3294]38   USE timing          ! Timing
[2715]39
[3]40   IMPLICIT NONE
41   PRIVATE
42
[508]43   PUBLIC   istate_init   ! routine called by step.F90
[3]44
45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
[4990]48   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
[888]49   !! $Id$
[2715]50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE istate_init
55      !!----------------------------------------------------------------------
56      !!                   ***  ROUTINE istate_init  ***
57      !!
[508]58      !! ** Purpose :   Initialization of the dynamics and tracer fields.
[3]59      !!----------------------------------------------------------------------
[7753]60      INTEGER ::   ji, jj, jk   ! dummy loop indices
[8568]61!!gm see comment further down
62      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
63!!gm end
[3294]64      !!----------------------------------------------------------------------
65      !
[8568]66      IF( ln_timing )   CALL timing_start('istate_init')
[3294]67      !
[5836]68      IF(lwp) WRITE(numout,*)
[7646]69      IF(lwp) WRITE(numout,*) 'istate_init : Initialization of the dynamics and tracers'
70      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[3]71
[7646]72!!gm  Why not include in the first call of dta_tsd ? 
73!!gm  probably associated with the use of internal damping...
[5836]74                     CALL dta_tsd_init        ! Initialisation of T & S input data
[7646]75!!gm to be moved in usrdef of C1D case
76!      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
77!!gm
[3]78
[7753]79      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
80      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
81      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
82      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
83
[15]84      IF( ln_rstart ) THEN                    ! Restart from a file
[3]85         !                                    ! -------------------
[6140]86         CALL rst_read                        ! Read the restart file
[7646]87         IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry
[6140]88         CALL day_init                        ! model calendar (using both namelist and restart infos)
[7646]89         !
90      ELSE                                    ! Start from rest
[3]91         !                                    ! ---------------
[7646]92         numror = 0                           ! define numror = 0 -> no restart file to read
93         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
94         CALL day_init                        ! model calendar (using both namelist and restart infos)
95         !                                    ! Initialization of ocean to zero
[508]96         !
[7646]97         IF( ln_tsd_init ) THEN               
98            CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
99            !
[7753]100            sshb(:,:)   = 0._wp               ! set the ocean at rest
101            ub  (:,:,:) = 0._wp
102            vb  (:,:,:) = 0._wp 
[7646]103            !
104         ELSE                                 ! user defined initial T and S
105            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
[3]106         ENDIF
[7753]107         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
108         sshn (:,:)     = sshb(:,:)   
109         un   (:,:,:)   = ub  (:,:,:)
110         vn   (:,:,:)   = vb  (:,:,:)
111         hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
[7646]112         CALL div_hor( 0 )                    ! compute interior hdivn value 
113!!gm                                    hdivn(:,:,:) = 0._wp
114
115!!gm POTENTIAL BUG :
116!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
117!!             as well as gdept and gdepw....   !!!!!
118!!      ===>>>>   probably a call to domvvl initialisation here....
119
120
[6140]121         !
[7646]122!!gm to be moved in usrdef of C1D case
123!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
[8568]124!            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
[7646]125!            CALL dta_uvd( nit000, zuvd )
126!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
127!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
[8568]128!            DEALLOCATE( zuvd )
[7646]129!         ENDIF
130         !
[6140]131!!gm This is to be changed !!!!
[7646]132!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
133!         IF( .NOT.ln_linssh ) THEN
134!            DO jk = 1, jpk
135!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
136!            END DO
137!         ENDIF
[6140]138!!gm
[2148]139         !
[7646]140      ENDIF 
[4354]141      !
142      ! Initialize "now" and "before" barotropic velocities:
[6140]143      ! Do it whatever the free surface method, these arrays being eventually used
[4354]144      !
[7753]145      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
146      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
[4354]147      !
[7646]148!!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
[4354]149      DO jk = 1, jpkm1
150         DO jj = 1, jpj
151            DO ji = 1, jpi
[6140]152               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
153               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
[4354]154               !
[6140]155               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
156               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
[4354]157            END DO
158         END DO
159      END DO
160      !
[7753]161      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
162      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
[4354]163      !
[7753]164      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
165      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
166      !
[8568]167      IF( ln_timing )   CALL timing_stop('istate_init')
[3294]168      !
[3]169   END SUBROUTINE istate_init
170
[7646]171   !!======================================================================
[3]172END MODULE istate
Note: See TracBrowser for help on using the repository browser.