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.

dtadyn.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90 @ 11027

Last change on this file since 11027 was 11027, checked in by acc, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Final renaming conversions and removal of temporary pointers. All non-AGRIF SETTE tests are passing (including test cases). AGRIF tests compile and link but segment on first call to Agrif_Regrid. NST changes are therefore a work in progress but nothing is broken that was not broken before
Property svn:keywords set to `Id`
File size: 42.1 KB

Rev	Line
[325]	1	MODULE dtadyn
	2	!!======================================================================
	3	!! * MODULE dtadyn *
[2528]	4	!! Off-line : interpolation of the physical fields
	5	!!======================================================================
	6	!! History : OPA ! 1992-01 (M. Imbard) Original code
	7	!! 8.0 ! 1998-04 (L.Bopp MA Foujols) slopes for isopyc.
	8	!! - ! 1998-05 (L. Bopp) read output of coupled run
	9	!! 8.2 ! 2001-01 (M. Levy et M. Benjelloul) add netcdf FORMAT
	10	!! NEMO 1.0 ! 2005-03 (O. Aumont and A. El Moussaoui) F90
	11	!! - ! 2005-12 (C. Ethe) Adapted for DEGINT
	12	!! 3.0 ! 2007-06 (C. Ethe) use of iom module
	13	!! 3.3 ! 2010-11 (C. Ethe) Full reorganization of the off-line: phasing with the on-line
[3294]	14	!! 3.4 ! 2011-05 (C. Ethe) Use of fldread
[2528]	15	!!----------------------------------------------------------------------
[325]	16
	17	!!----------------------------------------------------------------------
[3294]	18	!! dta_dyn_init : initialization, namelist read, and SAVEs control
[325]	19	!! dta_dyn : Interpolation of the fields
	20	!!----------------------------------------------------------------------
	21	USE oce ! ocean dynamics and tracers variables
[2528]	22	USE c1d ! 1D configuration: lk_c1d
	23	USE dom_oce ! ocean domain: variables
[7646]	24	USE domvvl ! variable volume
[2528]	25	USE zdf_oce ! ocean vertical physics: variables
	26	USE sbc_oce ! surface module: variables
[3294]	27	USE trc_oce ! share ocean/biogeo variables
[325]	28	USE phycst ! physical constants
[2528]	29	USE trabbl ! active tracer: bottom boundary layer
	30	USE ldfslp ! lateral diffusion: iso-neutral slopes
[7646]	31	USE sbcrnf ! river runoffs
	32	USE ldftra ! ocean tracer lateral physics
[2528]	33	USE zdfmxl ! vertical physics: mixed layer depth
	34	USE eosbn2 ! equation of state - Brunt Vaisala frequency
[325]	35	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[2528]	36	USE zpshde ! z-coord. with partial steps: horizontal derivatives
	37	USE in_out_manager ! I/O manager
	38	USE iom ! I/O library
[325]	39	USE lib_mpp ! distributed memory computing library
[3294]	40	USE prtctl ! print control
	41	USE fldread ! read input fields
	42	USE timing ! Timing
[7646]	43	USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc
[325]	44
	45	IMPLICIT NONE
	46	PRIVATE
	47
[10222]	48	PUBLIC dta_dyn_init ! called by opa.F90
	49	PUBLIC dta_dyn ! called by step.F90
	50	PUBLIC dta_dyn_sed_init ! called by opa.F90
	51	PUBLIC dta_dyn_sed ! called by step.F90
	52	PUBLIC dta_dyn_swp ! called by step.F90
[325]	53
[7646]	54	CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files
	55	LOGICAL :: ln_dynrnf !: read runoff data in file (T) or set to zero (F)
	56	LOGICAL :: ln_dynrnf_depth !: read runoff data in file (T) or set to zero (F)
	57	REAL(wp) :: fwbcorr
[325]	58
[7646]	59
	60	INTEGER , PARAMETER :: jpfld = 20 ! maximum number of fields to read
[3294]	61	INTEGER , SAVE :: jf_tem ! index of temperature
	62	INTEGER , SAVE :: jf_sal ! index of salinity
[7646]	63	INTEGER , SAVE :: jf_uwd ! index of u-transport
	64	INTEGER , SAVE :: jf_vwd ! index of v-transport
	65	INTEGER , SAVE :: jf_wwd ! index of v-transport
[3294]	66	INTEGER , SAVE :: jf_avt ! index of Kz
	67	INTEGER , SAVE :: jf_mld ! index of mixed layer deptht
	68	INTEGER , SAVE :: jf_emp ! index of water flux
[7646]	69	INTEGER , SAVE :: jf_empb ! index of water flux
[3294]	70	INTEGER , SAVE :: jf_qsr ! index of solar radiation
	71	INTEGER , SAVE :: jf_wnd ! index of wind speed
	72	INTEGER , SAVE :: jf_ice ! index of sea ice cover
[4570]	73	INTEGER , SAVE :: jf_rnf ! index of river runoff
[7646]	74	INTEGER , SAVE :: jf_fmf ! index of downward salt flux
[3294]	75	INTEGER , SAVE :: jf_ubl ! index of u-bbl coef
	76	INTEGER , SAVE :: jf_vbl ! index of v-bbl coef
[7646]	77	INTEGER , SAVE :: jf_div ! index of e3t
[325]	78
[7646]	79
	80	TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn ! structure of input fields (file informations, fields read)
[3294]	81	! !
	82	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta ! zonal isopycnal slopes
	83	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta ! meridional isopycnal slopes
	84	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta ! zonal diapycnal slopes
	85	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta ! meridional diapycnal slopes
[1735]	86
[7646]	87	INTEGER, SAVE :: nprevrec, nsecdyn
[325]	88
[343]	89	!!----------------------------------------------------------------------
[9598]	90	!! NEMO/OFF 4.0 , NEMO Consortium (2018)
[2528]	91	!! $Id$
[10068]	92	!! Software governed by the CeCILL license (see ./LICENSE)
[343]	93	!!----------------------------------------------------------------------
[325]	94	CONTAINS
	95
[10921]	96	SUBROUTINE dta_dyn( kt, Kbb, Kmm, Kaa )
[325]	97	!!----------------------------------------------------------------------
	98	!! * ROUTINE dta_dyn *
	99	!!
[3294]	100	!! ** Purpose : Prepares dynamics and physics fields from a NEMO run
	101	!! for an off-line simulation of passive tracers
[325]	102	!!
[3294]	103	!! ** Method : calculates the position of data
	104	!! - computes slopes if needed
	105	!! - interpolates data if needed
[2528]	106	!!----------------------------------------------------------------------
[10921]	107	INTEGER, INTENT(in) :: kt ! ocean time-step index
	108	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
[9212]	109	!
[7646]	110	INTEGER :: ji, jj, jk
[9212]	111	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zemp
	112	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdivtr
[325]	113	!!----------------------------------------------------------------------
[3294]	114	!
[9124]	115	IF( ln_timing ) CALL timing_start( 'dta_dyn')
[3294]	116	!
[7646]	117	nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
[3294]	118	!
[7646]	119	IF( kt == nit000 ) THEN ; nprevrec = 0
	120	ELSE ; nprevrec = sf_dyn(jf_tem)%nrec_a(2)
[325]	121	ENDIF
[7646]	122	CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==!
	123	!
[10922]	124	IF( l_ldfslp .AND. .NOT.lk_c1d ) CALL dta_dyn_slp( kt, Kbb, Kmm ) ! Computation of slopes
[7646]	125	!
[10921]	126	ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
	127	ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
[7646]	128	wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange
	129	fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+)
	130	fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction
	131	qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation
	132	emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P
	133	IF( ln_dynrnf ) THEN
	134	rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P
[10921]	135	IF( ln_dynrnf_depth .AND. .NOT. ln_linssh ) CALL dta_dyn_hrnf(Kmm)
[3294]	136	ENDIF
[325]	137	!
[10921]	138	uu(:,:,:,Kmm) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport
	139	vv(:,:,:,Kmm) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport
	140	ww(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport
[7646]	141	!
	142	IF( .NOT.ln_linssh ) THEN
[9212]	143	ALLOCATE( zemp(jpi,jpj) , zhdivtr(jpi,jpj,jpk) )
	144	zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:) ! effective u-transport
	145	emp_b (:,:) = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1) ! E-P
	146	zemp (:,:) = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr ) * tmask(:,:,1)
[10921]	147	CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa), e3t(:,:,:,Kaa) ) != ssh, vertical scale factor & vertical transport
[9212]	148	DEALLOCATE( zemp , zhdivtr )
[7646]	149	! Write in the tracer restart file
[9212]	150	! *********************************
[7646]	151	IF( lrst_trc ) THEN
[3827]	152	IF(lwp) WRITE(numout,*)
[9212]	153	IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp
	154	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[10921]	155	CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Kaa) )
	156	CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Kmm) )
[1501]	157	ENDIF
[3294]	158	ENDIF
[325]	159	!
[10921]	160	CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
[10954]	161	CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) ! now local thermal/haline expension ratio at T-points
	162	CALL bn2 ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) ! before Brunt-Vaisala frequency need for zdfmxl
[5131]	163
[10955]	164	rn2b(:,:,:) = rn2(:,:,:) ! needed for zdfmxl
	165	CALL zdf_mxl( kt, Kmm ) ! In any case, we need mxl
[2528]	166	!
[9019]	167	hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht
	168	avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient
[10213]	169	avs(:,:,:) = avt(:,:,:)
[7646]	170	!
[9019]	171	IF( ln_trabbl .AND. .NOT.lk_c1d ) THEN ! diffusive Bottom boundary layer param
	172	ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) ! bbl diffusive coef
	173	ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
	174	ENDIF
[2762]	175	!
[7646]	176	!
[10921]	177	CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
[7646]	178	!
[3294]	179	IF(ln_ctl) THEN ! print control
[10921]	180	CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask, kdim=jpk )
	181	CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn - : ', mask1=tmask, kdim=jpk )
	182	CALL prt_ctl(tab3d_1=uu(:,:,:,Kmm) , clinfo1=' uu(:,:,:,Kmm) - : ', mask1=umask, kdim=jpk )
	183	CALL prt_ctl(tab3d_1=vv(:,:,:,Kmm) , clinfo1=' vv(:,:,:,Kmm) - : ', mask1=vmask, kdim=jpk )
[11027]	184	CALL prt_ctl(tab3d_1=ww , clinfo1=' ww - : ', mask1=tmask, kdim=jpk )
[9440]	185	CALL prt_ctl(tab3d_1=avt , clinfo1=' kz - : ', mask1=tmask, kdim=jpk )
[7646]	186	CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
	187	CALL prt_ctl(tab3d_1=wslpi , clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
[2528]	188	ENDIF
	189	!
[9124]	190	IF( ln_timing ) CALL timing_stop( 'dta_dyn')
[3294]	191	!
[325]	192	END SUBROUTINE dta_dyn
	193
[2528]	194
[10921]	195	SUBROUTINE dta_dyn_init( Kbb, Kmm, Kaa )
[325]	196	!!----------------------------------------------------------------------
[3294]	197	!! * ROUTINE dta_dyn_init *
[325]	198	!!
[3294]	199	!! ** Purpose : Initialisation of the dynamical data
	200	!! ** Method : - read the data namdta_dyn namelist
[325]	201	!!----------------------------------------------------------------------
[10921]	202	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
	203	!
[3294]	204	INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code
	205	INTEGER :: ifpr ! dummy loop indice
	206	INTEGER :: jfld ! dummy loop arguments
	207	INTEGER :: inum, idv, idimv ! local integer
[4147]	208	INTEGER :: ios ! Local integer output status for namelist read
[7646]	209	INTEGER :: ji, jj, jk
	210	REAL(wp) :: zcoef
	211	INTEGER :: nkrnf_max
	212	REAL(wp) :: hrnf_max
[3294]	213	!!
[7646]	214	CHARACTER(len=100) :: cn_dir ! Root directory for location of core files
	215	TYPE(FLD_N), DIMENSION(jpfld) :: slf_d ! array of namelist informations on the fields to read
	216	TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp ! informations about the fields to be read
	217	TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt ! " "
	218	TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf ! " "
	219	TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf ! " "
	220	TYPE(FLD_N) :: sn_div ! informations about the fields to be read
[9212]	221	!!
[7646]	222	NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, &
[9212]	223	& sn_uwd, sn_vwd, sn_wwd, sn_emp, &
	224	& sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , &
	225	& sn_wnd, sn_ice, sn_fmf, &
	226	& sn_ubl, sn_vbl, sn_rnf, &
[7646]	227	& sn_empb, sn_div
[9212]	228	!!----------------------------------------------------------------------
[7646]	229	!
[4147]	230	REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
	231	READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
[9169]	232	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
[4147]	233	REWIND( numnam_cfg ) ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
	234	READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
[9169]	235	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
[4624]	236	IF(lwm) WRITE ( numond, namdta_dyn )
[3294]	237	! ! store namelist information in an array
	238	! ! Control print
[325]	239	IF(lwp) THEN
	240	WRITE(numout,*)
[3294]	241	WRITE(numout,*) 'dta_dyn : offline dynamics '
	242	WRITE(numout,*) '~~~~~~~ '
	243	WRITE(numout,*) ' Namelist namdta_dyn'
[7646]	244	WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf
	245	WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth
	246	WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr
[325]	247	WRITE(numout,*)
	248	ENDIF
[3294]	249	!
[7646]	250	jf_uwd = 1 ; jf_vwd = 2 ; jf_wwd = 3 ; jf_emp = 4 ; jf_avt = 5
	251	jf_tem = 6 ; jf_sal = 7 ; jf_mld = 8 ; jf_qsr = 9
	252	jf_wnd = 10 ; jf_ice = 11 ; jf_fmf = 12 ; jfld = jf_fmf
[3294]	253	!
[7646]	254	slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd
	255	slf_d(jf_emp) = sn_emp ; slf_d(jf_avt) = sn_avt
	256	slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld
	257	slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_ice) = sn_ice
	258	slf_d(jf_fmf) = sn_fmf
[3294]	259	!
[7646]	260	IF( .NOT.ln_linssh ) THEN
[9212]	261	jf_div = jfld + 1 ; jf_empb = jfld + 2 ; jfld = jf_empb
	262	slf_d(jf_div) = sn_div ; slf_d(jf_empb) = sn_empb
[7646]	263	ENDIF
	264	!
[9019]	265	IF( ln_trabbl ) THEN
[9212]	266	jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl
	267	slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl
[7646]	268	ENDIF
	269	!
[5385]	270	IF( ln_dynrnf ) THEN
[9212]	271	jf_rnf = jfld + 1 ; jfld = jf_rnf
	272	slf_d(jf_rnf) = sn_rnf
[4570]	273	ELSE
[7646]	274	rnf(:,:) = 0._wp
[4570]	275	ENDIF
	276
[3294]	277	ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure
[7646]	278	IF( ierr > 0 ) THEN
[3294]	279	CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN
	280	ENDIF
[5768]	281	! ! fill sf with slf_i and control print
	282	CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
[7646]	283	!
[3294]	284	! Open file for each variable to get his number of dimension
	285	DO ifpr = 1, jfld
[5768]	286	CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
	287	idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar
	288	idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar
	289	IF( sf_dyn(ifpr)%num /= 0 ) CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
	290	ierr1=0
[3294]	291	IF( idimv == 3 ) THEN ! 2D variable
	292	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 )
	293	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 )
	294	ELSE ! 3D variable
	295	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
	296	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
[2528]	297	ENDIF
[3294]	298	IF( ierr0 + ierr1 > 0 ) THEN
	299	CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN
	300	ENDIF
	301	END DO
[325]	302	!
[5836]	303	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! slopes
[3294]	304	IF( sf_dyn(jf_tem)%ln_tint ) THEN ! time interpolation
	305	ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2), &
	306	& wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
[7646]	307	!
	308	IF( ierr2 > 0 ) THEN
	309	CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN
	310	ENDIF
[3294]	311	ENDIF
[2528]	312	ENDIF
[7646]	313	!
	314	IF( .NOT.ln_linssh ) THEN
	315	IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND. & ! Restart: read in restart file
	316	iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
[10921]	317	IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation'
	318	CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Kmm) )
	319	CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Kbb) )
[7646]	320	ELSE
[10921]	321	IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation'
[7646]	322	CALL iom_open( 'restart', inum )
[10921]	323	CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Kmm) )
	324	CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Kbb) )
[7646]	325	CALL iom_close( inum ) ! close file
	326	ENDIF
	327	!
	328	DO jk = 1, jpkm1
[10921]	329	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
[7646]	330	ENDDO
[10921]	331	e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)
[7646]	332
	333	! Horizontal scale factor interpolations
	334	! --------------------------------------
[10921]	335	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
	336	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
[7646]	337
	338	! Vertical scale factor interpolations
	339	! ------------------------------------
[10921]	340	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' )
[7646]	341
[10921]	342	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
	343	e3u(:,:,:,Kbb) = e3u(:,:,:,Kmm)
	344	e3v(:,:,:,Kbb) = e3v(:,:,:,Kmm)
[7646]	345
	346	! t- and w- points depth
	347	! ----------------------
[10921]	348	gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
	349	gdepw(:,:,1,Kmm) = 0.0_wp
[7646]	350
	351	DO jk = 2, jpk
	352	DO jj = 1,jpj
	353	DO ji = 1,jpi
	354	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere
	355	! tmask = wmask, ie everywhere expect at jk = mikt
	356	! 1 for jk =
	357	! mikt
	358	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
[10921]	359	gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
	360	gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) &
	361	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
[7646]	362	END DO
	363	END DO
	364	END DO
	365
[10921]	366	gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
	367	gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
[7646]	368	!
[495]	369	ENDIF
[2715]	370	!
[7646]	371	IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed
	372	IF(lwp) WRITE(numout,*)
	373	IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
	374	CALL iom_open ( "runoffs", inum ) ! open file
	375	CALL iom_get ( inum, jpdom_data, 'rodepth', h_rnf ) ! read the river mouth array
	376	CALL iom_close( inum ) ! close file
	377	!
	378	nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied
	379	DO jj = 1, jpj
	380	DO ji = 1, jpi
	381	IF( h_rnf(ji,jj) > 0._wp ) THEN
	382	jk = 2
	383	DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
	384	END DO
	385	nk_rnf(ji,jj) = jk
	386	ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1
	387	ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj)
	388	ELSE
	389	CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' )
	390	WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
	391	ENDIF
	392	END DO
	393	END DO
	394	DO jj = 1, jpj ! set the associated depth
	395	DO ji = 1, jpi
	396	h_rnf(ji,jj) = 0._wp
	397	DO jk = 1, nk_rnf(ji,jj)
[10921]	398	h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
[7646]	399	END DO
	400	END DO
	401	END DO
	402	ELSE ! runoffs applied at the surface
	403	nk_rnf(:,:) = 1
[10921]	404	h_rnf (:,:) = e3t(:,:,1,Kmm)
[7646]	405	ENDIF
	406	nkrnf_max = MAXVAL( nk_rnf(:,:) )
	407	hrnf_max = MAXVAL( h_rnf(:,:) )
	408	IF( lk_mpp ) THEN
[10425]	409	CALL mpp_max( 'dtadyn', nkrnf_max ) ! max over the global domain
	410	CALL mpp_max( 'dtadyn', hrnf_max ) ! max over the global domain
[7646]	411	ENDIF
	412	IF(lwp) WRITE(numout,*) ' '
	413	IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,' max level : ', nkrnf_max
	414	IF(lwp) WRITE(numout,*) ' '
	415	!
[10921]	416	CALL dta_dyn( nit000, Kbb, Kmm, Kaa )
[2528]	417	!
[1501]	418	END SUBROUTINE dta_dyn_init
	419
[10921]	420	SUBROUTINE dta_dyn_sed( kt, Kmm )
[10222]	421	!!----------------------------------------------------------------------
	422	!! * ROUTINE dta_dyn *
	423	!!
	424	!! ** Purpose : Prepares dynamics and physics fields from a NEMO run
	425	!! for an off-line simulation of passive tracers
	426	!!
	427	!! ** Method : calculates the position of data
	428	!! - computes slopes if needed
	429	!! - interpolates data if needed
	430	!!----------------------------------------------------------------------
	431	INTEGER, INTENT(in) :: kt ! ocean time-step index
[10921]	432	INTEGER, INTENT(in) :: Kmm ! ocean time level index
[10222]	433	!
	434	!!----------------------------------------------------------------------
	435	!
	436	IF( ln_timing ) CALL timing_start( 'dta_dyn_sed')
	437	!
	438	nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
	439	!
	440	IF( kt == nit000 ) THEN ; nprevrec = 0
	441	ELSE ; nprevrec = sf_dyn(jf_tem)%nrec_a(2)
	442	ENDIF
	443	CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==!
	444	!
[10921]	445	ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
	446	ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
[10222]	447	!
[10921]	448	CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
[9212]	449
[10222]	450	IF(ln_ctl) THEN ! print control
[10921]	451	CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask, kdim=jpk )
	452	CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn - : ', mask1=tmask, kdim=jpk )
[10222]	453	ENDIF
	454	!
	455	IF( ln_timing ) CALL timing_stop( 'dta_dyn_sed')
	456	!
	457	END SUBROUTINE dta_dyn_sed
	458
	459
[10921]	460	SUBROUTINE dta_dyn_sed_init( Kmm )
[10222]	461	!!----------------------------------------------------------------------
	462	!! * ROUTINE dta_dyn_init *
	463	!!
	464	!! ** Purpose : Initialisation of the dynamical data
	465	!! ** Method : - read the data namdta_dyn namelist
	466	!!----------------------------------------------------------------------
[10921]	467	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
	468	!
[10222]	469	INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code
	470	INTEGER :: ifpr ! dummy loop indice
	471	INTEGER :: jfld ! dummy loop arguments
	472	INTEGER :: inum, idv, idimv ! local integer
	473	INTEGER :: ios ! Local integer output status for namelist read
	474	!!
	475	CHARACTER(len=100) :: cn_dir ! Root directory for location of core files
	476	TYPE(FLD_N), DIMENSION(2) :: slf_d ! array of namelist informations on the fields to read
	477	TYPE(FLD_N) :: sn_tem , sn_sal ! " "
	478	!!
	479	NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, &
	480	& sn_tem, sn_sal
	481	!!----------------------------------------------------------------------
	482	!
	483	REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
	484	READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
	485	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
	486	REWIND( numnam_cfg ) ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
	487	READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
	488	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
	489	IF(lwm) WRITE ( numond, namdta_dyn )
	490	! ! store namelist information in an array
	491	! ! Control print
	492	IF(lwp) THEN
	493	WRITE(numout,*)
	494	WRITE(numout,*) 'dta_dyn : offline dynamics '
	495	WRITE(numout,*) '~~~~~~~ '
	496	WRITE(numout,*) ' Namelist namdta_dyn'
	497	WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf
	498	WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth
	499	WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr
	500	WRITE(numout,*)
	501	ENDIF
	502	!
	503	jf_tem = 1 ; jf_sal = 2 ; jfld = jf_sal
	504	!
	505	slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal
	506	!
	507	ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure
	508	IF( ierr > 0 ) THEN
	509	CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN
	510	ENDIF
	511	! ! fill sf with slf_i and control print
	512	CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
	513	!
	514	! Open file for each variable to get his number of dimension
	515	DO ifpr = 1, jfld
	516	CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
	517	idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar
	518	idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar
	519	IF( sf_dyn(ifpr)%num /= 0 ) CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
	520	ierr1=0
	521	IF( idimv == 3 ) THEN ! 2D variable
	522	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 )
	523	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 )
	524	ELSE ! 3D variable
	525	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
	526	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
	527	ENDIF
	528	IF( ierr0 + ierr1 > 0 ) THEN
	529	CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN
	530	ENDIF
	531	END DO
	532	!
[10921]	533	CALL dta_dyn_sed( nit000, Kmm )
[10222]	534	!
	535	END SUBROUTINE dta_dyn_sed_init
	536
[10921]	537	SUBROUTINE dta_dyn_swp( kt, Kbb, Kmm, Kaa )
[7646]	538	!!---------------------------------------------------------------------
	539	!! * ROUTINE dta_dyn_swp *
	540	!!
[9212]	541	!! ** Purpose : Swap and the data and compute the vertical scale factor
	542	!! at U/V/W pointand the depht
[7646]	543	!!---------------------------------------------------------------------
[10921]	544	INTEGER, INTENT(in) :: kt ! time step
	545	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
[9212]	546	!
[7646]	547	INTEGER :: ji, jj, jk
	548	REAL(wp) :: zcoef
	549	!!---------------------------------------------------------------------
[6140]	550
[7646]	551	IF( kt == nit000 ) THEN
	552	IF(lwp) WRITE(numout,*)
	553	IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
	554	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	555	ENDIF
	556
[10921]	557	ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa)) ! before <-- now filtered
	558	ssh(:,:,Kmm) = ssh(:,:,Kaa)
[7646]	559
[10921]	560	e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)
[7646]	561
	562	! Reconstruction of all vertical scale factors at now and before time steps
	563	! =============================================================================
	564
	565	! Horizontal scale factor interpolations
	566	! --------------------------------------
[10921]	567	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
	568	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
[7646]	569
	570	! Vertical scale factor interpolations
	571	! ------------------------------------
[10921]	572	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' )
[7646]	573
[10921]	574	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
	575	e3u(:,:,:,Kbb) = e3u(:,:,:,Kmm)
	576	e3v(:,:,:,Kbb) = e3v(:,:,:,Kmm)
[7646]	577
	578	! t- and w- points depth
	579	! ----------------------
[10921]	580	gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
	581	gdepw(:,:,1,Kmm) = 0.0_wp
[9212]	582	!
[7646]	583	DO jk = 2, jpk
	584	DO jj = 1,jpj
	585	DO ji = 1,jpi
[9212]	586	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
[10921]	587	gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
	588	gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) &
	589	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
[9212]	590	END DO
	591	END DO
	592	END DO
	593	!
[10921]	594	gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
	595	gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
[7646]	596	!
	597	END SUBROUTINE dta_dyn_swp
[9212]	598
[7646]	599
	600	SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta )
[1501]	601	!!----------------------------------------------------------------------
[7646]	602	!! * ROUTINE dta_dyn_wzv *
	603	!!
[10921]	604	!! ** Purpose : compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity
[1501]	605	!!
[7646]	606	!! ** Method : Using the incompressibility hypothesis,
	607	!! - the ssh increment is computed by integrating the horizontal divergence
	608	!! and multiply by the time step.
[1501]	609	!!
[7646]	610	!! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly
	611	!! to the level thickness ( z-star case )
	612	!!
	613	!! - the vertical velocity is computed by integrating the horizontal divergence
	614	!! from the bottom to the surface minus the scale factor evolution.
	615	!! The boundary conditions are w=0 at the bottom (no flux)
	616	!!
[10921]	617	!! ** action : ssh(:,:,Kaa) / e3t(:,:,:,Kaa) / ww
[7646]	618	!!
	619	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[2528]	620	!!----------------------------------------------------------------------
[7646]	621	INTEGER, INTENT(in ) :: kt ! time-step
	622	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport
	623	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh
	624	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation
	625	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh
	626	REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor
[9212]	627	!
[7646]	628	INTEGER :: jk
	629	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
	630	REAL(wp) :: z2dt
	631	!!----------------------------------------------------------------------
[3294]	632	!
[7646]	633	z2dt = 2._wp * rdt
	634	!
	635	zhdiv(:,:) = 0._wp
	636	DO jk = 1, jpkm1
	637	zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk)
	638	END DO
	639	! ! Sea surface elevation time-stepping
	640	pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
	641	! !
	642	! ! After acale factors at t-points ( z_star coordinate )
	643	DO jk = 1, jpkm1
	644	pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
	645	END DO
	646	!
	647	END SUBROUTINE dta_dyn_ssh
	648
	649
[10921]	650	SUBROUTINE dta_dyn_hrnf( Kmm )
[7646]	651	!!----------------------------------------------------------------------
	652	!! * ROUTINE sbc_rnf *
[1501]	653	!!
[7646]	654	!! ** Purpose : update the horizontal divergence with the runoff inflow
	655	!!
	656	!! ** Method :
	657	!! CAUTION : rnf is positive (inflow) decreasing the
	658	!! divergence and expressed in m/s
	659	!!
	660	!! ** Action : phdivn decreased by the runoff inflow
[2528]	661	!!----------------------------------------------------------------------
[7646]	662	!!
[10921]	663	INTEGER, INTENT(in) :: Kmm ! ocean time level index
	664	!
	665	INTEGER :: ji, jj, jk ! dummy loop indices
[7646]	666	!!----------------------------------------------------------------------
[2528]	667	!
[7646]	668	DO jj = 1, jpj ! update the depth over which runoffs are distributed
	669	DO ji = 1, jpi
	670	h_rnf(ji,jj) = 0._wp
	671	DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres
[10921]	672	h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) ! to the bottom of the relevant grid box
[1501]	673	END DO
[7646]	674	END DO
[2528]	675	END DO
[5836]	676	!
[7646]	677	END SUBROUTINE dta_dyn_hrnf
	678
	679
	680
[10922]	681	SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm )
[7646]	682	!!---------------------------------------------------------------------
	683	!! * ROUTINE dta_dyn_slp *
	684	!!
	685	!! ** Purpose : Computation of slope
	686	!!
	687	!!---------------------------------------------------------------------
	688	INTEGER, INTENT(in) :: kt ! time step
[10922]	689	INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices
[7646]	690	!
	691	INTEGER :: ji, jj ! dummy loop indices
	692	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
	693	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
	694	INTEGER :: iswap
[9212]	695	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zuslp, zvslp, zwslpi, zwslpj
	696	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts
[7646]	697	!!---------------------------------------------------------------------
	698	!
	699	IF( sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace)
	700	IF( kt == nit000 ) THEN
	701	IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
	702	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature
	703	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity
	704	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef.
[10922]	705	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]	706	uslpdta (:,:,:,1) = zuslp (:,:,:)
	707	vslpdta (:,:,:,1) = zvslp (:,:,:)
	708	wslpidta(:,:,:,1) = zwslpi(:,:,:)
	709	wslpjdta(:,:,:,1) = zwslpj(:,:,:)
	710	!
	711	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature
	712	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity
	713	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef.
[10922]	714	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]	715	uslpdta (:,:,:,2) = zuslp (:,:,:)
	716	vslpdta (:,:,:,2) = zvslp (:,:,:)
	717	wslpidta(:,:,:,2) = zwslpi(:,:,:)
	718	wslpjdta(:,:,:,2) = zwslpj(:,:,:)
	719	ELSE
	720	!
	721	iswap = 0
	722	IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 ) iswap = 1
	723	IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 ) THEN ! read/update the after data
	724	IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
	725	uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data
	726	vslpdta (:,:,:,1) = vslpdta (:,:,:,2)
	727	wslpidta(:,:,:,1) = wslpidta(:,:,:,2)
	728	wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2)
	729	!
	730	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature
	731	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity
	732	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef.
[10922]	733	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]	734	!
	735	uslpdta (:,:,:,2) = zuslp (:,:,:)
	736	vslpdta (:,:,:,2) = zvslp (:,:,:)
	737	wslpidta(:,:,:,2) = zwslpi(:,:,:)
	738	wslpjdta(:,:,:,2) = zwslpj(:,:,:)
	739	ENDIF
	740	ENDIF
	741	ENDIF
	742	!
	743	IF( sf_dyn(jf_tem)%ln_tint ) THEN
	744	ztinta = REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp ) &
	745	& / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
	746	ztintb = 1. - ztinta
	747	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
	748	uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2)
	749	vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2)
	750	wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2)
	751	wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2)
	752	ENDIF
	753	ELSE
	754	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
	755	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
	756	avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coef.
[10922]	757	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]	758	!
	759	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
	760	uslp (:,:,:) = zuslp (:,:,:)
	761	vslp (:,:,:) = zvslp (:,:,:)
	762	wslpi(:,:,:) = zwslpi(:,:,:)
	763	wslpj(:,:,:) = zwslpj(:,:,:)
	764	ENDIF
	765	ENDIF
	766	!
	767	END SUBROUTINE dta_dyn_slp
[1501]	768
[9212]	769
[10922]	770	SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm )
[1501]	771	!!---------------------------------------------------------------------
[3294]	772	!! * ROUTINE dta_dyn_slp *
[1501]	773	!!
[9212]	774	!! ** Purpose : Computation of slope
[1501]	775	!!---------------------------------------------------------------------
[3294]	776	INTEGER , INTENT(in ) :: kt ! time step
	777	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! temperature/salinity
	778	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: puslp ! zonal isopycnal slopes
	779	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pvslp ! meridional isopycnal slopes
	780	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpi ! zonal diapycnal slopes
	781	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpj ! meridional diapycnal slopes
[10922]	782	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
[1501]	783	!!---------------------------------------------------------------------
[9212]	784	!
[7646]	785	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
[5836]	786	CALL eos ( pts, rhd, rhop, gdept_0(:,:,:) )
[10954]	787	CALL eos_rab( pts, rab_n, Kmm ) ! now local thermal/haline expension ratio at T-points
	788	CALL bn2 ( pts, rab_n, rn2, Kmm ) ! now Brunt-Vaisala
[5131]	789
[6140]	790	! Partial steps: before Horizontal DErivative
	791	IF( ln_zps .AND. .NOT. ln_isfcav) &
[10954]	792	& CALL zps_hde ( kt, Kmm, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient
[6140]	793	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
	794	IF( ln_zps .AND. ln_isfcav) &
[10954]	795	& CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
[7646]	796	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
[4990]	797
[10955]	798	rn2b(:,:,:) = rn2(:,:,:) ! needed for zdfmxl
	799	CALL zdf_mxl( kt, Kmm ) ! mixed layer depth
[10922]	800	CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm ) ! slopes
[7646]	801	puslp (:,:,:) = uslp (:,:,:)
	802	pvslp (:,:,:) = vslp (:,:,:)
	803	pwslpi(:,:,:) = wslpi(:,:,:)
	804	pwslpj(:,:,:) = wslpj(:,:,:)
[5836]	805	ELSE
	806	puslp (:,:,:) = 0. ! to avoid warning when compiling
	807	pvslp (:,:,:) = 0.
	808	pwslpi(:,:,:) = 0.
	809	pwslpj(:,:,:) = 0.
	810	ENDIF
[2528]	811	!
[7646]	812	END SUBROUTINE compute_slopes
[9212]	813
[2528]	814	!!======================================================================
[325]	815	END MODULE dtadyn

Note: See TracBrowser for help on using the repository browser.

Download in other formats: