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 @ 13207

Last change on this file since 13207 was 12489, checked in by davestorkey, 4 years ago

Preparation for new timestepping scheme #2390.
Main changes:

Initial euler timestep now handled in stp and not in TRA/DYN routines.
Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
Renaming of a few miscellaneous parameters, eg. atfp -> rn_atfp (namelist parameter used everywhere) and rau0 -> rho0.

This version gives bit-comparable results to the previous version of the trunk.

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
[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
	67	INTEGER , SAVE :: jf_wwd ! index of v-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
	124	ELSE ; nprevrec = sf_dyn(jf_tem)%nrec_a(2)
[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' )
[7646]	285	!
[3294]	286	! Open file for each variable to get his number of dimension
	287	DO ifpr = 1, jfld
[12377]	288	CALL fld_def( sf_dyn(ifpr) )
	289	CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num )
[5768]	290	idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar
	291	idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar
[12377]	292	CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
[5768]	293	ierr1=0
[3294]	294	IF( idimv == 3 ) THEN ! 2D variable
	295	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 )
	296	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 )
	297	ELSE ! 3D variable
	298	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
	299	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
[2528]	300	ENDIF
[3294]	301	IF( ierr0 + ierr1 > 0 ) THEN
	302	CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN
	303	ENDIF
	304	END DO
[325]	305	!
[5836]	306	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! slopes
[3294]	307	IF( sf_dyn(jf_tem)%ln_tint ) THEN ! time interpolation
	308	ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2), &
	309	& wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
[7646]	310	!
	311	IF( ierr2 > 0 ) THEN
	312	CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN
	313	ENDIF
[3294]	314	ENDIF
[2528]	315	ENDIF
[7646]	316	!
	317	IF( .NOT.ln_linssh ) THEN
	318	IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND. & ! Restart: read in restart file
	319	iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
[12377]	320	IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation'
	321	CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Kmm) )
	322	CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Kbb) )
[7646]	323	ELSE
[12377]	324	IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation'
[7646]	325	CALL iom_open( 'restart', inum )
[12377]	326	CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Kmm) )
	327	CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Kbb) )
[7646]	328	CALL iom_close( inum ) ! close file
	329	ENDIF
	330	!
	331	DO jk = 1, jpkm1
[12377]	332	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
[7646]	333	ENDDO
[12377]	334	e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)
[7646]	335
	336	! Horizontal scale factor interpolations
	337	! --------------------------------------
[12377]	338	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
	339	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
[7646]	340
	341	! Vertical scale factor interpolations
	342	! ------------------------------------
[12377]	343	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' )
[7646]	344
[12377]	345	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
	346	e3u(:,:,:,Kbb) = e3u(:,:,:,Kmm)
	347	e3v(:,:,:,Kbb) = e3v(:,:,:,Kmm)
[7646]	348
	349	! t- and w- points depth
	350	! ----------------------
[12377]	351	gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
	352	gdepw(:,:,1,Kmm) = 0.0_wp
[7646]	353
[12377]	354	DO_3D_11_11( 2, jpk )
	355	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere
	356	! tmask = wmask, ie everywhere expect at jk = mikt
	357	! 1 for jk =
	358	! mikt
	359	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
	360	gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
	361	gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) &
	362	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
	363	END_3D
[7646]	364
[12377]	365	gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
	366	gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
[7646]	367	!
[495]	368	ENDIF
[2715]	369	!
[7646]	370	IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed
	371	IF(lwp) WRITE(numout,*)
	372	IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
	373	CALL iom_open ( "runoffs", inum ) ! open file
	374	CALL iom_get ( inum, jpdom_data, 'rodepth', h_rnf ) ! read the river mouth array
	375	CALL iom_close( inum ) ! close file
	376	!
	377	nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied
[12377]	378	DO_2D_11_11
	379	IF( h_rnf(ji,jj) > 0._wp ) THEN
	380	jk = 2
	381	DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
[7646]	382	END DO
[12377]	383	nk_rnf(ji,jj) = jk
	384	ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1
	385	ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj)
	386	ELSE
	387	CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' )
	388	WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
	389	ENDIF
	390	END_2D
	391	DO_2D_11_11
	392	h_rnf(ji,jj) = 0._wp
	393	DO jk = 1, nk_rnf(ji,jj)
	394	h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
[7646]	395	END DO
[12377]	396	END_2D
[7646]	397	ELSE ! runoffs applied at the surface
	398	nk_rnf(:,:) = 1
[12377]	399	h_rnf (:,:) = e3t(:,:,1,Kmm)
[7646]	400	ENDIF
	401	nkrnf_max = MAXVAL( nk_rnf(:,:) )
	402	hrnf_max = MAXVAL( h_rnf(:,:) )
	403	IF( lk_mpp ) THEN
[10425]	404	CALL mpp_max( 'dtadyn', nkrnf_max ) ! max over the global domain
	405	CALL mpp_max( 'dtadyn', hrnf_max ) ! max over the global domain
[7646]	406	ENDIF
	407	IF(lwp) WRITE(numout,*) ' '
	408	IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,' max level : ', nkrnf_max
	409	IF(lwp) WRITE(numout,*) ' '
	410	!
[12377]	411	CALL dta_dyn( nit000, Kbb, Kmm, Kaa )
[2528]	412	!
[1501]	413	END SUBROUTINE dta_dyn_init
	414
[12377]	415	SUBROUTINE dta_dyn_sed( kt, Kmm )
[10222]	416	!!----------------------------------------------------------------------
	417	!! * ROUTINE dta_dyn *
	418	!!
	419	!! ** Purpose : Prepares dynamics and physics fields from a NEMO run
	420	!! for an off-line simulation of passive tracers
	421	!!
	422	!! ** Method : calculates the position of data
	423	!! - computes slopes if needed
	424	!! - interpolates data if needed
	425	!!----------------------------------------------------------------------
	426	INTEGER, INTENT(in) :: kt ! ocean time-step index
[12377]	427	INTEGER, INTENT(in) :: Kmm ! ocean time level index
[10222]	428	!
	429	!!----------------------------------------------------------------------
	430	!
	431	IF( ln_timing ) CALL timing_start( 'dta_dyn_sed')
	432	!
	433	nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
	434	!
	435	IF( kt == nit000 ) THEN ; nprevrec = 0
	436	ELSE ; nprevrec = sf_dyn(jf_tem)%nrec_a(2)
	437	ENDIF
	438	CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==!
	439	!
[12377]	440	ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
	441	ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
[10222]	442	!
[12377]	443	CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
[9212]	444
[12377]	445	IF(sn_cfctl%l_prtctl) THEN ! print control
	446	CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask, kdim=jpk )
	447	CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn - : ', mask1=tmask, kdim=jpk )
[10222]	448	ENDIF
	449	!
	450	IF( ln_timing ) CALL timing_stop( 'dta_dyn_sed')
	451	!
	452	END SUBROUTINE dta_dyn_sed
	453
	454
[12377]	455	SUBROUTINE dta_dyn_sed_init( Kmm )
[10222]	456	!!----------------------------------------------------------------------
	457	!! * ROUTINE dta_dyn_init *
	458	!!
	459	!! ** Purpose : Initialisation of the dynamical data
	460	!! ** Method : - read the data namdta_dyn namelist
	461	!!----------------------------------------------------------------------
[12377]	462	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
	463	!
[10222]	464	INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code
	465	INTEGER :: ifpr ! dummy loop indice
	466	INTEGER :: jfld ! dummy loop arguments
	467	INTEGER :: inum, idv, idimv ! local integer
	468	INTEGER :: ios ! Local integer output status for namelist read
	469	!!
	470	CHARACTER(len=100) :: cn_dir ! Root directory for location of core files
	471	TYPE(FLD_N), DIMENSION(2) :: slf_d ! array of namelist informations on the fields to read
	472	TYPE(FLD_N) :: sn_tem , sn_sal ! " "
	473	!!
	474	NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, &
	475	& sn_tem, sn_sal
	476	!!----------------------------------------------------------------------
	477	!
	478	READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
[11536]	479	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' )
[10222]	480	READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
[11536]	481	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' )
[10222]	482	IF(lwm) WRITE ( numond, namdta_dyn )
	483	! ! store namelist information in an array
	484	! ! Control print
	485	IF(lwp) THEN
	486	WRITE(numout,*)
	487	WRITE(numout,*) 'dta_dyn : offline dynamics '
	488	WRITE(numout,*) '~~~~~~~ '
	489	WRITE(numout,*) ' Namelist namdta_dyn'
	490	WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf
	491	WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth
	492	WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr
	493	WRITE(numout,*)
	494	ENDIF
	495	!
	496	jf_tem = 1 ; jf_sal = 2 ; jfld = jf_sal
	497	!
	498	slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal
	499	!
	500	ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure
	501	IF( ierr > 0 ) THEN
	502	CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN
	503	ENDIF
	504	! ! fill sf with slf_i and control print
	505	CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
	506	!
	507	! Open file for each variable to get his number of dimension
	508	DO ifpr = 1, jfld
[12377]	509	CALL fld_def( sf_dyn(ifpr) )
	510	CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num )
[10222]	511	idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar
	512	idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar
[12377]	513	CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
[10222]	514	ierr1=0
	515	IF( idimv == 3 ) THEN ! 2D variable
	516	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 )
	517	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 )
	518	ELSE ! 3D variable
	519	ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 )
	520	IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
	521	ENDIF
	522	IF( ierr0 + ierr1 > 0 ) THEN
	523	CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN
	524	ENDIF
	525	END DO
	526	!
[12377]	527	CALL dta_dyn_sed( nit000, Kmm )
[10222]	528	!
	529	END SUBROUTINE dta_dyn_sed_init
	530
[12377]	531	SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa )
[7646]	532	!!---------------------------------------------------------------------
	533	!! * ROUTINE dta_dyn_swp *
	534	!!
[12377]	535	!! ** Purpose : Asselin time filter of now SSH
[7646]	536	!!---------------------------------------------------------------------
[12377]	537	INTEGER, INTENT(in) :: kt ! time step
	538	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
[9212]	539	!
[7646]	540	!!---------------------------------------------------------------------
[6140]	541
[7646]	542	IF( kt == nit000 ) THEN
	543	IF(lwp) WRITE(numout,*)
[12377]	544	IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height'
	545	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
[7646]	546	ENDIF
	547
[12489]	548	ssh(:,:,Kmm) = ssh(:,:,Kmm) + rn_atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa))
[7646]	549
[12377]	550	!! Do we also need to time filter e3t??
	551	!
	552	END SUBROUTINE dta_dyn_atf
	553
	554	SUBROUTINE dta_dyn_sf_interp( kt, Kmm )
	555	!!---------------------------------------------------------------------
	556	!! * ROUTINE dta_dyn_sf_interp *
	557	!!
	558	!! ** Purpose : Calculate scale factors at U/V/W points and depths
	559	!! given the after e3t field
	560	!!---------------------------------------------------------------------
	561	INTEGER, INTENT(in) :: kt ! time step
	562	INTEGER, INTENT(in) :: Kmm ! ocean time level indices
	563	!
	564	INTEGER :: ji, jj, jk
	565	REAL(wp) :: zcoef
	566	!!---------------------------------------------------------------------
[7646]	567
	568	! Horizontal scale factor interpolations
	569	! --------------------------------------
[12377]	570	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
	571	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
[7646]	572
	573	! Vertical scale factor interpolations
	574	! ------------------------------------
[12377]	575	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' )
[7646]	576
	577	! t- and w- points depth
	578	! ----------------------
[12377]	579	gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
	580	gdepw(:,:,1,Kmm) = 0.0_wp
[9212]	581	!
[12377]	582	DO_3D_11_11( 2, jpk )
	583	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
	584	gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
	585	gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) &
	586	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
	587	END_3D
[9212]	588	!
[12377]	589	END SUBROUTINE dta_dyn_sf_interp
[7646]	590
	591	SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta )
[1501]	592	!!----------------------------------------------------------------------
[7646]	593	!! * ROUTINE dta_dyn_wzv *
	594	!!
[12377]	595	!! ** Purpose : compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity
[1501]	596	!!
[7646]	597	!! ** Method : Using the incompressibility hypothesis,
	598	!! - the ssh increment is computed by integrating the horizontal divergence
	599	!! and multiply by the time step.
[1501]	600	!!
[7646]	601	!! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly
	602	!! to the level thickness ( z-star case )
	603	!!
	604	!! - the vertical velocity is computed by integrating the horizontal divergence
	605	!! from the bottom to the surface minus the scale factor evolution.
	606	!! The boundary conditions are w=0 at the bottom (no flux)
	607	!!
[12377]	608	!! ** action : ssh(:,:,Kaa) / e3t(:,:,:,Kaa) / ww
[7646]	609	!!
	610	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[2528]	611	!!----------------------------------------------------------------------
[7646]	612	INTEGER, INTENT(in ) :: kt ! time-step
	613	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport
	614	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh
	615	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation
	616	REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh
	617	REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor
[9212]	618	!
[7646]	619	INTEGER :: jk
	620	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
	621	REAL(wp) :: z2dt
	622	!!----------------------------------------------------------------------
[3294]	623	!
[12489]	624	z2dt = 2._wp * rn_Dt
[7646]	625	!
	626	zhdiv(:,:) = 0._wp
	627	DO jk = 1, jpkm1
	628	zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk)
	629	END DO
	630	! ! Sea surface elevation time-stepping
[12489]	631	pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rho0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
[7646]	632	! !
	633	! ! After acale factors at t-points ( z_star coordinate )
	634	DO jk = 1, jpkm1
	635	pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
	636	END DO
	637	!
	638	END SUBROUTINE dta_dyn_ssh
	639
	640
[12377]	641	SUBROUTINE dta_dyn_hrnf( Kmm )
[7646]	642	!!----------------------------------------------------------------------
	643	!! * ROUTINE sbc_rnf *
[1501]	644	!!
[7646]	645	!! ** Purpose : update the horizontal divergence with the runoff inflow
	646	!!
	647	!! ** Method :
	648	!! CAUTION : rnf is positive (inflow) decreasing the
	649	!! divergence and expressed in m/s
	650	!!
	651	!! ** Action : phdivn decreased by the runoff inflow
[2528]	652	!!----------------------------------------------------------------------
[7646]	653	!!
[12377]	654	INTEGER, INTENT(in) :: Kmm ! ocean time level index
	655	!
	656	INTEGER :: ji, jj, jk ! dummy loop indices
[7646]	657	!!----------------------------------------------------------------------
[2528]	658	!
[12377]	659	DO_2D_11_11
	660	h_rnf(ji,jj) = 0._wp
	661	DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres
	662	h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) ! to the bottom of the relevant grid box
	663	END DO
	664	END_2D
[5836]	665	!
[7646]	666	END SUBROUTINE dta_dyn_hrnf
	667
	668
	669
[12377]	670	SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm )
[7646]	671	!!---------------------------------------------------------------------
	672	!! * ROUTINE dta_dyn_slp *
	673	!!
	674	!! ** Purpose : Computation of slope
	675	!!
	676	!!---------------------------------------------------------------------
	677	INTEGER, INTENT(in) :: kt ! time step
[12377]	678	INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices
[7646]	679	!
	680	INTEGER :: ji, jj ! dummy loop indices
	681	REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation
	682	REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation
	683	INTEGER :: iswap
[9212]	684	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zuslp, zvslp, zwslpi, zwslpj
	685	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts
[7646]	686	!!---------------------------------------------------------------------
	687	!
	688	IF( sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace)
	689	IF( kt == nit000 ) THEN
	690	IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
	691	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature
	692	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity
	693	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef.
[12377]	694	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]	695	uslpdta (:,:,:,1) = zuslp (:,:,:)
	696	vslpdta (:,:,:,1) = zvslp (:,:,:)
	697	wslpidta(:,:,:,1) = zwslpi(:,:,:)
	698	wslpjdta(:,:,:,1) = zwslpj(:,:,:)
	699	!
	700	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature
	701	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity
	702	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef.
[12377]	703	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]	704	uslpdta (:,:,:,2) = zuslp (:,:,:)
	705	vslpdta (:,:,:,2) = zvslp (:,:,:)
	706	wslpidta(:,:,:,2) = zwslpi(:,:,:)
	707	wslpjdta(:,:,:,2) = zwslpj(:,:,:)
	708	ELSE
	709	!
	710	iswap = 0
	711	IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 ) iswap = 1
	712	IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 ) THEN ! read/update the after data
	713	IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
	714	uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data
	715	vslpdta (:,:,:,1) = vslpdta (:,:,:,2)
	716	wslpidta(:,:,:,1) = wslpidta(:,:,:,2)
	717	wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2)
	718	!
	719	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature
	720	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity
	721	avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef.
[12377]	722	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]	723	!
	724	uslpdta (:,:,:,2) = zuslp (:,:,:)
	725	vslpdta (:,:,:,2) = zvslp (:,:,:)
	726	wslpidta(:,:,:,2) = zwslpi(:,:,:)
	727	wslpjdta(:,:,:,2) = zwslpj(:,:,:)
	728	ENDIF
	729	ENDIF
	730	ENDIF
	731	!
	732	IF( sf_dyn(jf_tem)%ln_tint ) THEN
	733	ztinta = REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp ) &
	734	& / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
	735	ztintb = 1. - ztinta
	736	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
	737	uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2)
	738	vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2)
	739	wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2)
	740	wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2)
	741	ENDIF
	742	ELSE
	743	zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature
	744	zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity
	745	avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coef.
[12377]	746	CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm )
[7646]	747	!
	748	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
	749	uslp (:,:,:) = zuslp (:,:,:)
	750	vslp (:,:,:) = zvslp (:,:,:)
	751	wslpi(:,:,:) = zwslpi(:,:,:)
	752	wslpj(:,:,:) = zwslpj(:,:,:)
	753	ENDIF
	754	ENDIF
	755	!
	756	END SUBROUTINE dta_dyn_slp
[1501]	757
[9212]	758
[12377]	759	SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm )
[1501]	760	!!---------------------------------------------------------------------
[3294]	761	!! * ROUTINE dta_dyn_slp *
[1501]	762	!!
[9212]	763	!! ** Purpose : Computation of slope
[1501]	764	!!---------------------------------------------------------------------
[3294]	765	INTEGER , INTENT(in ) :: kt ! time step
	766	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! temperature/salinity
	767	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: puslp ! zonal isopycnal slopes
	768	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pvslp ! meridional isopycnal slopes
	769	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpi ! zonal diapycnal slopes
	770	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpj ! meridional diapycnal slopes
[12377]	771	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
[1501]	772	!!---------------------------------------------------------------------
[9212]	773	!
[7646]	774	IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace)
[5836]	775	CALL eos ( pts, rhd, rhop, gdept_0(:,:,:) )
[12377]	776	CALL eos_rab( pts, rab_n, Kmm ) ! now local thermal/haline expension ratio at T-points
	777	CALL bn2 ( pts, rab_n, rn2, Kmm ) ! now Brunt-Vaisala
[5131]	778
[6140]	779	! Partial steps: before Horizontal DErivative
	780	IF( ln_zps .AND. .NOT. ln_isfcav) &
[12377]	781	& CALL zps_hde ( kt, Kmm, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient
[6140]	782	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
	783	IF( ln_zps .AND. ln_isfcav) &
[12377]	784	& CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
[7646]	785	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
[4990]	786
[12377]	787	rn2b(:,:,:) = rn2(:,:,:) ! needed for zdfmxl
	788	CALL zdf_mxl( kt, Kmm ) ! mixed layer depth
	789	CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm ) ! slopes
[7646]	790	puslp (:,:,:) = uslp (:,:,:)
	791	pvslp (:,:,:) = vslp (:,:,:)
	792	pwslpi(:,:,:) = wslpi(:,:,:)
	793	pwslpj(:,:,:) = wslpj(:,:,:)
[5836]	794	ELSE
	795	puslp (:,:,:) = 0. ! to avoid warning when compiling
	796	pvslp (:,:,:) = 0.
	797	pwslpi(:,:,:) = 0.
	798	pwslpj(:,:,:) = 0.
	799	ENDIF
[2528]	800	!
[7646]	801	END SUBROUTINE compute_slopes
[9212]	802
[2528]	803	!!======================================================================
[325]	804	END MODULE dtadyn

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

Download in other formats: