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/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OFF – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OFF/dtadyn.F90 @ 12960

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

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

Download in other formats: