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

Last change on this file since 10213 was 10213, checked in by aumont, 6 years ago
corrects bugs described in Tickets 2144 and 2145
Property svn:keywords set to `Id`
File size: 35.4 KB

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

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

Download in other formats: