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_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OFF – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13747_ENHANCE-04_dford_OBSOP_BGC/src/OFF/dtadyn.F90 @ 14049

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

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

Download in other formats: