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/trunk/src/OFF – NEMO

Context Navigation

source: NEMO/trunk/src/OFF/dtadyn.F90 @ 14255

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

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

Download in other formats: