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.

istate.F90 in branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 6418

Last change on this file since 6418 was 6020, checked in by timgraham, 9 years ago
Merged with head of trunk
Property svn:keywords set to `Id`
File size: 24.4 KB

Rev	Line
[3]	1	MODULE istate
	2	!!======================================================================
	3	!! * MODULE istate *
	4	!! Ocean state : initial state setting
	5	!!=====================================================================
[2104]	6	!! History : OPA ! 1989-12 (P. Andrich) Original code
	7	!! 5.0 ! 1991-11 (G. Madec) rewritting
	8	!! 6.0 ! 1996-01 (G. Madec) terrain following coordinates
	9	!! 8.0 ! 2001-09 (M. Levy, M. Ben Jelloul) istate_eel
	10	!! 8.0 ! 2001-09 (M. Levy, M. Ben Jelloul) istate_uvg
	11	!! NEMO 1.0 ! 2003-08 (G. Madec, C. Talandier) F90: Free form, modules + EEL R5
	12	!! - ! 2004-05 (A. Koch-Larrouy) istate_gyre
	13	!! 2.0 ! 2006-07 (S. Masson) distributed restart using iom
	14	!! 3.3 ! 2010-10 (C. Ethe) merge TRC-TRA
[3294]	15	!! 3.4 ! 2011-04 (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn
[508]	16	!!----------------------------------------------------------------------
[3]	17
	18	!!----------------------------------------------------------------------
	19	!! istate_init : initial state setting
	20	!! istate_tem : analytical profile for initial Temperature
	21	!! istate_sal : analytical profile for initial Salinity
	22	!! istate_eel : initial state setting of EEL R5 configuration
[93]	23	!! istate_gyre : initial state setting of GYRE configuration
[3]	24	!! istate_uvg : initial velocity in geostropic balance
	25	!!----------------------------------------------------------------------
	26	USE oce ! ocean dynamics and active tracers
	27	USE dom_oce ! ocean space and time domain
[4245]	28	USE c1d ! 1D vertical configuration
[2104]	29	USE daymod ! calendar
	30	USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine)
[6020]	31	USE ldftra ! lateral physics: ocean active tracers
[3]	32	USE zdf_oce ! ocean vertical physics
	33	USE phycst ! physical constants
[3294]	34	USE dtatsd ! data temperature and salinity (dta_tsd routine)
[4245]	35	USE dtauvd ! data: U & V current (dta_uvd routine)
[593]	36	USE domvvl ! varying vertical mesh
[4990]	37	!
	38	USE in_out_manager ! I/O manager
	39	USE iom ! I/O library
[2715]	40	USE lib_mpp ! MPP library
[3680]	41	USE restart ! restart
[3294]	42	USE wrk_nemo ! Memory allocation
	43	USE timing ! Timing
[2715]	44
[3]	45	IMPLICIT NONE
	46	PRIVATE
	47
[508]	48	PUBLIC istate_init ! routine called by step.F90
[3]	49
	50	!! * Substitutions
	51	# include "domzgr_substitute.h90"
	52	# include "vectopt_loop_substitute.h90"
	53	!!----------------------------------------------------------------------
[4990]	54	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
[888]	55	!! $Id$
[2715]	56	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	57	!!----------------------------------------------------------------------
	58	CONTAINS
	59
	60	SUBROUTINE istate_init
	61	!!----------------------------------------------------------------------
	62	!! * ROUTINE istate_init *
	63	!!
[508]	64	!! ** Purpose : Initialization of the dynamics and tracer fields.
[3]	65	!!----------------------------------------------------------------------
[5163]	66	INTEGER :: ji, jj, jk ! dummy loop indices
	67	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zuvd ! U & V data workspace
[3294]	68	!!----------------------------------------------------------------------
	69	!
[4990]	70	IF( nn_timing == 1 ) CALL timing_start('istate_init')
[3294]	71	!
[3]	72
[6020]	73	IF(lwp) WRITE(numout,*)
[508]	74	IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
	75	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[3]	76
[6020]	77	CALL dta_tsd_init ! Initialisation of T & S input data
	78	IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data
[3]	79
[5163]	80	rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk
	81	rn2b (:,:,: ) = 0._wp ; rn2 (:,:,: ) = 0._wp ! set one for all to 0 at levels 1 and jpk
	82	tsa (:,:,:,:) = 0._wp ! set one for all to 0 at level jpk
	83	rab_b(:,:,:,:) = 0._wp ; rab_n(:,:,:,:) = 0._wp ! set one for all to 0 at level jpk
[3294]	84
[15]	85	IF( ln_rstart ) THEN ! Restart from a file
[3]	86	! ! -------------------
	87	CALL rst_read ! Read the restart file
[1130]	88	CALL day_init ! model calendar (using both namelist and restart infos)
[3]	89	ELSE
	90	! ! Start from rest
	91	! ! ---------------
[1130]	92	numror = 0 ! define numror = 0 -> no restart file to read
[3]	93	neuler = 0 ! Set time-step indicator at nit000 (euler forward)
[1130]	94	CALL day_init ! model calendar (using both namelist and restart infos)
[508]	95	! ! Initialization of ocean to zero
[3294]	96	! before fields ! now fields
	97	sshb (:,:) = 0._wp ; sshn (:,:) = 0._wp
	98	ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp
	99	vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp
[6020]	100	hdivn(:,:,:) = 0._wp
[508]	101	!
[3]	102	IF( cp_cfg == 'eel' ) THEN
[2104]	103	CALL istate_eel ! EEL configuration : start from pre-defined U,V T-S fields
[434]	104	ELSEIF( cp_cfg == 'gyre' ) THEN
[2104]	105	CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields
[4245]	106	ELSE ! Initial T-S, U-V fields read in files
	107	IF ( ln_tsd_init ) THEN ! read 3D T and S data at nit000
	108	CALL dta_tsd( nit000, tsb )
	109	tsn(:,:,:,:) = tsb(:,:,:,:)
	110	!
	111	ELSE ! Initial T-S fields defined analytically
	112	CALL istate_t_s
	113	ENDIF
	114	IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
[6020]	115	CALL wrk_alloc( jpi,jpj,jpk,2, zuvd )
[4245]	116	CALL dta_uvd( nit000, zuvd )
	117	ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:)
	118	vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:)
[6020]	119	CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd )
[4245]	120	ENDIF
[3]	121	ENDIF
[2148]	122	!
	123	! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
	124	IF( lk_vvl ) THEN
	125	DO jk = 1, jpk
	126	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
[6020]	127	END DO
[2148]	128	ENDIF
	129	!
[3]	130	ENDIF
[1438]	131	!
[4354]	132	!
	133	! Initialize "now" and "before" barotropic velocities:
	134	! Do it whatever the free surface method, these arrays
	135	! being eventually used
	136	!
	137	!
[6020]	138	un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
	139	ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
[4354]	140	!
	141	DO jk = 1, jpkm1
	142	DO jj = 1, jpj
	143	DO ji = 1, jpi
	144	un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
	145	vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
	146	!
[4370]	147	ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
	148	vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
[4354]	149	END DO
	150	END DO
	151	END DO
	152	!
[4370]	153	un_b(:,:) = un_b(:,:) * hur (:,:)
	154	vn_b(:,:) = vn_b(:,:) * hvr (:,:)
[4354]	155	!
[4370]	156	ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
	157	vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
[4354]	158	!
	159	!
[4990]	160	IF( nn_timing == 1 ) CALL timing_stop('istate_init')
[3294]	161	!
[3]	162	END SUBROUTINE istate_init
	163
[4990]	164
[3294]	165	SUBROUTINE istate_t_s
[3]	166	!!---------------------------------------------------------------------
[3294]	167	!! * ROUTINE istate_t_s *
[3]	168	!!
	169	!! ** Purpose : Intialization of the temperature field with an
	170	!! analytical profile or a file (i.e. in EEL configuration)
	171	!!
[3294]	172	!! ** Method : - temperature: use Philander analytic profile
	173	!! - salinity : use to a constant value 35.5
[3]	174	!!
	175	!! References : Philander ???
	176	!!----------------------------------------------------------------------
[6020]	177	INTEGER :: ji, jj, jk
	178	REAL(wp) :: zsal = 35.50_wp
[3]	179	!!----------------------------------------------------------------------
[508]	180	!
[3]	181	IF(lwp) WRITE(numout,*)
[3294]	182	IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
	183	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)'
[2104]	184	!
[3]	185	DO jk = 1, jpk
[3294]	186	tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) ) &
	187	& + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.) ) * tmask(:,:,jk)
	188	tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
[3]	189	END DO
[3294]	190	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
	191	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
[2104]	192	!
[3294]	193	END SUBROUTINE istate_t_s
[3]	194
[4990]	195
[3]	196	SUBROUTINE istate_eel
	197	!!----------------------------------------------------------------------
	198	!! * ROUTINE istate_eel *
	199	!!
	200	!! ** Purpose : Initialization of the dynamics and tracers for EEL R5
	201	!! configuration (channel with or without a topographic bump)
	202	!!
	203	!! ** Method : - set temprature field
	204	!! - set salinity field
	205	!! - set velocity field including horizontal divergence
	206	!! and relative vorticity fields
	207	!!----------------------------------------------------------------------
[6020]	208	USE divhor ! hor. divergence (div_hor routine)
[473]	209	USE iom
[4990]	210	!
[3]	211	INTEGER :: inum ! temporary logical unit
	212	INTEGER :: ji, jj, jk ! dummy loop indices
[479]	213	INTEGER :: ijloc
[508]	214	REAL(wp) :: zh1, zh2, zslope, zcst, zfcor ! temporary scalars
[2104]	215	REAL(wp) :: zt1 = 15._wp ! surface temperature value (EEL R5)
	216	REAL(wp) :: zt2 = 5._wp ! bottom temperature value (EEL R5)
	217	REAL(wp) :: zsal = 35.0_wp ! constant salinity (EEL R2, R5 and R6)
	218	REAL(wp) :: zueel = 0.1_wp ! constant uniform zonal velocity (EEL R5)
[508]	219	REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain
[3]	220	!!----------------------------------------------------------------------
[4990]	221	!
[3]	222	SELECT CASE ( jp_cfg )
	223	! ! ====================
	224	CASE ( 5 ) ! EEL R5 configuration
	225	! ! ====================
[2104]	226	!
[3]	227	! set temperature field with a linear profile
	228	! -------------------------------------------
	229	IF(lwp) WRITE(numout,*)
	230	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
	231	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[2104]	232	!
[4292]	233	zh1 = gdept_1d( 1 )
	234	zh2 = gdept_1d(jpkm1)
[2104]	235	!
[3]	236	zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
	237	zcst = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
[2104]	238	!
[3]	239	DO jk = 1, jpk
[3294]	240	tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
	241	tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
[3]	242	END DO
[2104]	243	!
[3294]	244	IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &
	245	& 1 , jpi , 5 , 1 , jpk , &
	246	& 1 , 1. , numout )
[2104]	247	!
[3]	248	! set salinity field to a constant value
	249	! --------------------------------------
	250	IF(lwp) WRITE(numout,*)
	251	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
	252	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[2104]	253	!
[3294]	254	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
	255	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
[2104]	256	!
[6020]	257	! set the dynamics: U,V, hdiv (and ssh if necessary)
[3]	258	! ----------------
	259	! Start EEL5 configuration with barotropic geostrophic velocities
	260	! according the sshb and sshn SSH imposed.
[479]	261	! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
	262	! we use the Coriolis frequency at mid-channel.
	263	ub(:,:,:) = zueel * umask(:,:,:)
[3]	264	un(:,:,:) = ub(:,:,:)
[479]	265	ijloc = mj0(INT(jpjglo-1)/2)
	266	zfcor = ff(1,ijloc)
[2104]	267	!
[3]	268	DO jj = 1, jpjglo
[479]	269	zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)zueele1t(1,1)*zfcor/grav
[3]	270	END DO
[2104]	271	!
[479]	272	IF(lwp) THEN
	273	WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
	274	WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
	275	WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
	276	ENDIF
[2104]	277	!
[3]	278	DO jj = 1, nlcj
	279	DO ji = 1, nlci
	280	sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
	281	END DO
	282	END DO
	283	sshb(nlci+1:jpi, : ) = 0.e0 ! set to zero extra mpp columns
	284	sshb( : ,nlcj+1:jpj) = 0.e0 ! set to zero extra mpp rows
[2104]	285	!
[3]	286	sshn(:,:) = sshb(:,:) ! set now ssh to the before value
[2104]	287	!
[593]	288	IF( nn_rstssh /= 0 ) THEN
[2104]	289	nn_rstssh = 0 ! hand-made initilization of ssh
[593]	290	CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
[558]	291	ENDIF
[2104]	292	!
[6020]	293	!!gm Check here call to div_hor should not be necessary
	294	!!gm div_hor call runoffs not sure they are defined at that level
	295	CALL div_hor( nit000 ) ! horizontal divergence and relative vorticity (curl)
[3]	296	! N.B. the vertical velocity will be computed from the horizontal divergence field
	297	! in istate by a call to wzv routine
	298
	299
	300	! ! ==========================
	301	CASE ( 2 , 6 ) ! EEL R2 or R6 configuration
	302	! ! ==========================
[2104]	303	!
[3]	304	! set temperature field with a NetCDF file
	305	! ----------------------------------------
	306	IF(lwp) WRITE(numout,*)
	307	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
	308	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[2104]	309	!
[473]	310	CALL iom_open ( 'eel.initemp', inum )
[3294]	311	CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
[473]	312	CALL iom_close( inum )
[2104]	313	!
[3294]	314	tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb
[2104]	315	!
[3294]	316	IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &
	317	& 1 , jpi , 5 , 1 , jpk , &
	318	& 1 , 1. , numout )
[2104]	319	!
[3]	320	! set salinity field to a constant value
	321	! --------------------------------------
	322	IF(lwp) WRITE(numout,*)
	323	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
	324	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[2104]	325	!
[3294]	326	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
	327	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
[2104]	328	!
[3]	329	! ! ===========================
	330	CASE DEFAULT ! NONE existing configuration
	331	! ! ===========================
[473]	332	WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
	333	CALL ctl_stop( ctmp1 )
[2104]	334	!
[3]	335	END SELECT
[2104]	336	!
[3]	337	END SUBROUTINE istate_eel
	338
	339
[93]	340	SUBROUTINE istate_gyre
	341	!!----------------------------------------------------------------------
	342	!! * ROUTINE istate_gyre *
	343	!!
	344	!! ** Purpose : Initialization of the dynamics and tracers for GYRE
	345	!! configuration (double gyre with rotated domain)
	346	!!
	347	!! ** Method : - set temprature field
[6020]	348	!! - set salinity field
[93]	349	!!----------------------------------------------------------------------
[473]	350	INTEGER :: ji, jj, jk ! dummy loop indices
[508]	351	INTEGER :: inum ! temporary logical unit
	352	INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization
[93]	353	!!----------------------------------------------------------------------
[4990]	354	!
[434]	355	SELECT CASE ( ntsinit)
[4990]	356	!
[434]	357	CASE ( 0 ) ! analytical T/S profil deduced from LEVITUS
	358	IF(lwp) WRITE(numout,*)
	359	IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
	360	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[4990]	361	!
[434]	362	DO jk = 1, jpk
	363	DO jj = 1, jpj
	364	DO ji = 1, jpi
[3294]	365	tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 ) ) &
[434]	366	& * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2 &
	367	& + ( 15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) ) &
	368	& - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.) &
	369	& + 7. * (1500. - fsdept(ji,jj,jk)) / 1500. ) &
	370	& * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
[3294]	371	tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
	372	tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
[434]	373
[3294]	374	tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 ) ) &
[434]	375	& * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2 &
	376	& + ( 35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000. &
	377	& - 1.62 * TANH( (fsdept(ji,jj,jk) - 60. ) / 650. ) &
	378	& + 0.2 * TANH( (fsdept(ji,jj,jk) - 35. ) / 100. ) &
	379	& + 0.2 * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.) ) &
	380	& * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
[3294]	381	tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
	382	tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
[434]	383	END DO
[93]	384	END DO
	385	END DO
[4990]	386	!
[434]	387	CASE ( 1 ) ! T/S data fields read in dta_tem.nc/data_sal.nc files
	388	IF(lwp) WRITE(numout,*)
	389	IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
	390	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	391	IF(lwp) WRITE(numout,*) ' NetCDF FORMAT'
	392
	393	! Read temperature field
	394	! ----------------------
[473]	395	CALL iom_open ( 'data_tem', inum )
[3294]	396	CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) )
[473]	397	CALL iom_close( inum )
[434]	398
[3294]	399	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
	400	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
[434]	401
	402	! Read salinity field
	403	! -------------------
[473]	404	CALL iom_open ( 'data_sal', inum )
[3294]	405	CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) )
[473]	406	CALL iom_close( inum )
[434]	407
[3294]	408	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
	409	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
[4990]	410	!
[434]	411	END SELECT
[4990]	412	!
[93]	413	IF(lwp) THEN
	414	WRITE(numout,*)
	415	WRITE(numout,*) ' Initial temperature and salinity profiles:'
[4292]	416	WRITE(numout, "(9x,' level gdept_1d temperature salinity ')" )
	417	WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )
[93]	418	ENDIF
[4990]	419	!
[93]	420	END SUBROUTINE istate_gyre
	421
[4990]	422
[3]	423	SUBROUTINE istate_uvg
	424	!!----------------------------------------------------------------------
	425	!! * ROUTINE istate_uvg *
	426	!!
	427	!! ** Purpose : Compute the geostrophic velocities from (tn,sn) fields
	428	!!
	429	!! ** Method : Using the hydrostatic hypothesis the now hydrostatic
	430	!! pressure is computed by integrating the in-situ density from the
	431	!! surface to the bottom.
	432	!! p=integral [ rau*g dz ]
	433	!!----------------------------------------------------------------------
[6020]	434	USE divhor ! hor. divergence (div_hor routine)
[3]	435	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[4990]	436	!
[3]	437	INTEGER :: ji, jj, jk ! dummy loop indices
[508]	438	REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars
[3294]	439	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
[3]	440	!!----------------------------------------------------------------------
[3294]	441	!
[6020]	442	CALL wrk_alloc( jpi,jpj,jpk, zprn)
[3294]	443	!
[3]	444	IF(lwp) WRITE(numout,*)
	445	IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
	446	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	447
	448	! Compute the now hydrostatic pressure
	449	! ------------------------------------
	450
[15]	451	zalfg = 0.5 * grav * rau0
[508]	452
	453	zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value
[3]	454
[508]	455	DO jk = 2, jpkm1 ! Vertical integration from the surface
[3]	456	zprn(:,:,jk) = zprn(:,:,jk-1) &
[359]	457	& + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
[3]	458	END DO
	459
	460	! Compute geostrophic balance
	461	! ---------------------------
	462	DO jk = 1, jpkm1
	463	DO jj = 2, jpjm1
	464	DO ji = fs_2, fs_jpim1 ! vertor opt.
	465	zmsv = 1. / MAX( umask(ji-1,jj+1,jk) + umask(ji ,jj+1,jk) &
	466	+ umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) , 1. )
	467	zphv = ( zprn(ji ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1) &
	468	+ ( zprn(ji+1,jj+1,jk) - zprn(ji ,jj+1,jk) ) * umask(ji ,jj+1,jk) / e1u(ji ,jj+1) &
	469	+ ( zprn(ji ,jj ,jk) - zprn(ji-1,jj ,jk) ) * umask(ji-1,jj ,jk) / e1u(ji-1,jj ) &
	470	+ ( zprn(ji+1,jj ,jk) - zprn(ji ,jj ,jk) ) * umask(ji ,jj ,jk) / e1u(ji ,jj )
	471	zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
	472
	473	zmsu = 1. / MAX( vmask(ji+1,jj ,jk) + vmask(ji ,jj ,jk) &
	474	+ vmask(ji+1,jj-1,jk) + vmask(ji ,jj-1,jk) , 1. )
	475	zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj ,jk) ) * vmask(ji+1,jj ,jk) / e2v(ji+1,jj ) &
	476	+ ( zprn(ji ,jj+1,jk) - zprn(ji ,jj ,jk) ) * vmask(ji ,jj ,jk) / e2v(ji ,jj ) &
	477	+ ( zprn(ji+1,jj ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1) &
	478	+ ( zprn(ji ,jj ,jk) - zprn(ji ,jj-1,jk) ) * vmask(ji ,jj-1,jk) / e2v(ji ,jj-1)
	479	zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
	480
	481	! Compute the geostrophic velocities
	482	un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji ,jj-1) )
	483	vn(ji,jj,jk) = 2. * zphv / ( ff(ji,jj) + ff(ji-1,jj ) )
	484	END DO
	485	END DO
	486	END DO
	487
	488	IF(lwp) WRITE(numout,*) ' we force to zero bottom velocity'
	489
	490	! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
	491	! to have a zero bottom velocity
	492
	493	DO jk = 1, jpkm1
	494	un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
	495	vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
	496	END DO
	497
	498	CALL lbc_lnk( un, 'U', -1. )
	499	CALL lbc_lnk( vn, 'V', -1. )
	500
	501	ub(:,:,:) = un(:,:,:)
	502	vb(:,:,:) = vn(:,:,:)
	503
[508]	504	!
[6020]	505	!!gm Check here call to div_hor should not be necessary
	506	!!gm div_hor call runoffs not sure they are defined at that level
	507	CALL div_hor( nit000 ) ! now horizontal divergence
[2715]	508	!
[6020]	509	CALL wrk_dealloc( jpi,jpj,jpk, zprn)
	510	!
[3]	511	END SUBROUTINE istate_uvg
	512
	513	!!=====================================================================
	514	END MODULE istate

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

Download in other formats: