source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/istate.F90 @ 13219

Last change on this file since 13219 was 13219, checked in by smasson, 4 months ago

better e3: update with trunk@13218 see #2385

  • Property svn:keywords set to Id
File size: 9.2 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 dtatsd         ! data temperature and salinity   (dta_tsd routine)
27   USE dtauvd         ! data: U & V current             (dta_uvd routine)
[593]28   USE domvvl          ! varying vertical mesh
[7646]29   USE wet_dry         ! wetting and drying (needed for wad_istate)
30   USE usrdef_istate   ! User defined initial state
[4990]31   !
32   USE in_out_manager  ! I/O manager
33   USE iom             ! I/O library
[2715]34   USE lib_mpp         ! MPP library
[3680]35   USE restart         ! restart
[13219]36#if defined key_agrif
37   USE agrif_oce_interp
38   USE agrif_oce
39#endif   
[2715]40
[3]41   IMPLICIT NONE
42   PRIVATE
43
[508]44   PUBLIC   istate_init   ! routine called by step.F90
[3]45
46   !! * Substitutions
[12377]47#  include "do_loop_substitute.h90"
[12680]48#  include "domzgr_substitute.h90"
[3]49   !!----------------------------------------------------------------------
[9598]50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]51   !! $Id$
[10068]52   !! Software governed by the CeCILL license (see ./LICENSE)
[3]53   !!----------------------------------------------------------------------
54CONTAINS
55
[12377]56   SUBROUTINE istate_init( Kbb, Kmm, Kaa )
[3]57      !!----------------------------------------------------------------------
58      !!                   ***  ROUTINE istate_init  ***
59      !!
[508]60      !! ** Purpose :   Initialization of the dynamics and tracer fields.
[3]61      !!----------------------------------------------------------------------
[12377]62      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices
63      !
[7753]64      INTEGER ::   ji, jj, jk   ! dummy loop indices
[12732]65      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table  !!st patch to use gdept subtitute
[9019]66!!gm see comment further down
67      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
68!!gm end
[3294]69      !!----------------------------------------------------------------------
70      !
[5836]71      IF(lwp) WRITE(numout,*)
[7646]72      IF(lwp) WRITE(numout,*) 'istate_init : Initialization of the dynamics and tracers'
73      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[3]74
[7646]75!!gm  Why not include in the first call of dta_tsd ? 
76!!gm  probably associated with the use of internal damping...
[13219]77       CALL dta_tsd_init        ! Initialisation of T & S input data
[7646]78!!gm to be moved in usrdef of C1D case
79!      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
80!!gm
[3]81
[7753]82      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
83      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
[12377]84      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
[7753]85      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
[9453]86#if defined key_agrif
[12377]87      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization
88      vv   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization   
[9453]89#endif
[7753]90
[13219]91#if defined key_agrif
92      IF ( (.NOT.Agrif_root()).AND.ln_init_chfrpar ) THEN
93         numror = 0                           ! define numror = 0 -> no restart file to read
94         ln_1st_euler = .true.                ! Set time-step indicator at nit000 (euler forward)
95         CALL day_init 
96         CALL agrif_istate( Kbb, Kmm, Kaa )   ! Interp from parent
97         !
98         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) 
99         ssh (:,:,Kmm)     = ssh(:,:,Kbb)
100         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
101         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
102      ELSE
103#endif
[15]104      IF( ln_rstart ) THEN                    ! Restart from a file
[3]105         !                                    ! -------------------
[12377]106         CALL rst_read( Kbb, Kmm )            ! Read the restart file
[6140]107         CALL day_init                        ! model calendar (using both namelist and restart infos)
[7646]108         !
109      ELSE                                    ! Start from rest
[3]110         !                                    ! ---------------
[7646]111         numror = 0                           ! define numror = 0 -> no restart file to read
[12724]112         l_1st_euler = .true.                 ! Set time-step indicator at nit000 (euler forward)
[7646]113         CALL day_init                        ! model calendar (using both namelist and restart infos)
114         !                                    ! Initialization of ocean to zero
[508]115         !
[7646]116         IF( ln_tsd_init ) THEN               
[12377]117            CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000
[7646]118            !
[12377]119            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest
[13219]120            uu  (:,:,:,Kbb) = 0._wp
121            vv  (:,:,:,Kbb) = 0._wp 
122            !
[10499]123            IF( ll_wd ) THEN
[12377]124               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD
[10499]125               !
126               ! Apply minimum wetdepth criterion
127               !
[12377]128               DO_2D_11_11
129                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
130                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
131                  ENDIF
132               END_2D
[10499]133            ENDIF 
[13219]134             !
[7646]135         ELSE                                 ! user defined initial T and S
[12911]136            DO jk = 1, jpk
[12731]137               zgdept(:,:,jk) = gdept(:,:,jk,Kbb)
138            END DO
139            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
[3]140         ENDIF
[12377]141         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
142         ssh (:,:,Kmm)     = ssh(:,:,Kbb)   
143         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
144         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
[7646]145
146!!gm POTENTIAL BUG :
[12377]147!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
[12731]148!!             as well as gdept_ and gdepw_....   !!!!!
[7646]149!!      ===>>>>   probably a call to domvvl initialisation here....
150
151
[6140]152         !
[7646]153!!gm to be moved in usrdef of C1D case
154!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
[9019]155!            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
[7646]156!            CALL dta_uvd( nit000, zuvd )
[12377]157!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
158!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
[9019]159!            DEALLOCATE( zuvd )
[7646]160!         ENDIF
161         !
[6140]162!!gm This is to be changed !!!!
[12377]163!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here
[7646]164!         IF( .NOT.ln_linssh ) THEN
165!            DO jk = 1, jpk
[12377]166!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
[7646]167!            END DO
168!         ENDIF
[6140]169!!gm
[2148]170         !
[7646]171      ENDIF 
[13219]172#if defined key_agrif
173      ENDIF
174#endif
[4354]175      !
176      ! Initialize "now" and "before" barotropic velocities:
[6140]177      ! Do it whatever the free surface method, these arrays being eventually used
[4354]178      !
[12377]179      uu_b(:,:,Kmm) = 0._wp   ;   vv_b(:,:,Kmm) = 0._wp
180      uu_b(:,:,Kbb) = 0._wp   ;   vv_b(:,:,Kbb) = 0._wp
[4354]181      !
[12377]182!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
183      DO_3D_11_11( 1, jpkm1 )
184         uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
185         vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
186         !
187         uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
188         vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
189      END_3D
[4354]190      !
[12377]191      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
192      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
[4354]193      !
[12377]194      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb)
195      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb)
[7753]196      !
[3]197   END SUBROUTINE istate_init
198
[7646]199   !!======================================================================
[3]200END MODULE istate
Note: See TracBrowser for help on using the repository browser.