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/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 3229

Last change on this file since 3229 was 3186, checked in by smasson, 13 years ago
dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo
Property svn:keywords set to `Id`
File size: 25.2 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
[2977]	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
[2104]	28	USE daymod ! calendar
	29	USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine)
[3]	30	USE ldftra_oce ! ocean active tracers: lateral physics
	31	USE zdf_oce ! ocean vertical physics
	32	USE phycst ! physical constants
[2977]	33	USE dtatsd ! data temperature and salinity (dta_tsd routine)
[3]	34	USE restart ! ocean restart (rst_read routine)
[508]	35	USE in_out_manager ! I/O manager
[2104]	36	USE iom ! I/O library
[544]	37	USE zpshde ! partial step: hor. derivative (zps_hde routine)
	38	USE eosbn2 ! equation of state (eos bn2 routine)
[593]	39	USE domvvl ! varying vertical mesh
	40	USE dynspg_oce ! pressure gradient schemes
	41	USE dynspg_flt ! pressure gradient schemes
	42	USE dynspg_exp ! pressure gradient schemes
	43	USE dynspg_ts ! pressure gradient schemes
[2715]	44	USE lib_mpp ! MPP library
[3186]	45	USE wrk_nemo ! Memory allocation
[3162]	46	USE timing ! Timing
[2715]	47
[3]	48	IMPLICIT NONE
	49	PRIVATE
	50
[508]	51	PUBLIC istate_init ! routine called by step.F90
[3]	52
	53	!! * Substitutions
	54	# include "domzgr_substitute.h90"
	55	# include "vectopt_loop_substitute.h90"
	56	!!----------------------------------------------------------------------
[2287]	57	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	58	!! $Id$
[2715]	59	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	60	!!----------------------------------------------------------------------
	61	CONTAINS
	62
	63	SUBROUTINE istate_init
	64	!!----------------------------------------------------------------------
	65	!! * ROUTINE istate_init *
	66	!!
[508]	67	!! ** Purpose : Initialization of the dynamics and tracer fields.
[3]	68	!!----------------------------------------------------------------------
[2148]	69	! - ML - needed for initialization of e3t_b
	70	INTEGER :: jk ! dummy loop indice
[3162]	71	!!----------------------------------------------------------------------
	72	!
	73	IF( nn_timing == 1 ) CALL timing_start('istate_init')
	74	!
[3]	75
[508]	76	IF(lwp) WRITE(numout,*)
	77	IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
	78	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[3]	79
[2977]	80	CALL dta_tsd_init ! Initialisation of T & S input data
[3]	81
[2977]	82	rhd (:,:,: ) = 0.e0
	83	rhop (:,:,: ) = 0.e0
	84	rn2 (:,:,: ) = 0.e0
	85	tsa (:,:,:,:) = 0.e0
	86
[15]	87	IF( ln_rstart ) THEN ! Restart from a file
[3]	88	! ! -------------------
	89	neuler = 1 ! Set time-step indicator at nit000 (leap-frog)
	90	CALL rst_read ! Read the restart file
[3116]	91	! ! define e3u_b, e3v_b from e3t_b read in restart file
	92	CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
[1130]	93	CALL day_init ! model calendar (using both namelist and restart infos)
[3]	94	ELSE
	95	! ! Start from rest
	96	! ! ---------------
[1130]	97	numror = 0 ! define numror = 0 -> no restart file to read
[3]	98	neuler = 0 ! Set time-step indicator at nit000 (euler forward)
[1130]	99	CALL day_init ! model calendar (using both namelist and restart infos)
[508]	100	! ! Initialization of ocean to zero
[3116]	101	! before fields ! now fields
	102	sshb (:,:) = 0._wp ; sshn (:,:) = 0._wp
	103	ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp
	104	vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp
	105	rotb (:,:,:) = 0._wp ; rotn (:,:,:) = 0._wp
	106	hdivb(:,:,:) = 0._wp ; hdivn(:,:,:) = 0._wp
[508]	107	!
[3116]	108	! ! define e3u_b, e3v_b from e3t_b initialized in domzgr
	109	CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
	110	!
	111	IF( cp_cfg == 'eel' ) THEN
[2104]	112	CALL istate_eel ! EEL configuration : start from pre-defined U,V T-S fields
[434]	113	ELSEIF( cp_cfg == 'gyre' ) THEN
[2104]	114	CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields
[2977]	115	ELSEIF( ln_tsd_init ) THEN ! Initial T-S fields read in files
	116	CALL dta_tsd( nit000, tsb ) ! read 3D T and S data at nit000
	117	tsn(:,:,:,:) = tsb(:,:,:,:)
	118	!
	119	ELSE ! Initial T-S fields defined analytically
	120	CALL istate_t_s
[3]	121	ENDIF
[2104]	122	!
[2082]	123	CALL eos( tsb, rhd, rhop ) ! before potential and in situ densities
[2236]	124	#if ! defined key_c1d
	125	IF( ln_zps ) CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv, & ! zps: before hor. gradient
	126	& rhd, gru , grv ) ! of t,s,rd at ocean bottom
	127	#endif
[2148]	128	!
	129	! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
	130	IF( lk_vvl ) THEN
	131	DO jk = 1, jpk
	132	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
	133	ENDDO
	134	ENDIF
	135	!
[3]	136	ENDIF
[1438]	137	!
[2104]	138	IF( lk_agrif ) THEN ! read free surface arrays in restart file
[593]	139	IF( ln_rstart ) THEN
	140	IF( lk_dynspg_flt ) CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields
[1438]	141	! ! gcx, gcxb for agrif_opa_init
	142	ENDIF ! explicit case not coded yet with AGRIF
[1200]	143	ENDIF
[508]	144	!
[3162]	145	IF( nn_timing == 1 ) CALL timing_stop('istate_init')
	146	!
[3]	147	END SUBROUTINE istate_init
	148
[2977]	149	SUBROUTINE istate_t_s
[3]	150	!!---------------------------------------------------------------------
[2977]	151	!! * ROUTINE istate_t_s *
[3]	152	!!
	153	!! ** Purpose : Intialization of the temperature field with an
	154	!! analytical profile or a file (i.e. in EEL configuration)
	155	!!
[2977]	156	!! ** Method : - temperature: use Philander analytic profile
	157	!! - salinity : use to a constant value 35.5
[3]	158	!!
	159	!! References : Philander ???
	160	!!----------------------------------------------------------------------
[2977]	161	INTEGER :: ji, jj, jk
	162	REAL(wp) :: zsal = 35.50
[3]	163	!!----------------------------------------------------------------------
[508]	164	!
[3]	165	IF(lwp) WRITE(numout,*)
[2977]	166	IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
	167	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)'
[2104]	168	!
[3]	169	DO jk = 1, jpk
[2977]	170	tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) ) &
	171	& + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.) ) * tmask(:,:,jk)
	172	tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
[3]	173	END DO
[2977]	174	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
	175	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
[2104]	176	!
[2977]	177	END SUBROUTINE istate_t_s
[3]	178
	179
	180	SUBROUTINE istate_eel
	181	!!----------------------------------------------------------------------
	182	!! * ROUTINE istate_eel *
	183	!!
	184	!! ** Purpose : Initialization of the dynamics and tracers for EEL R5
	185	!! configuration (channel with or without a topographic bump)
	186	!!
	187	!! ** Method : - set temprature field
	188	!! - set salinity field
	189	!! - set velocity field including horizontal divergence
	190	!! and relative vorticity fields
	191	!!----------------------------------------------------------------------
	192	USE divcur ! hor. divergence & rel. vorticity (div_cur routine)
[473]	193	USE iom
[3]	194
	195	INTEGER :: inum ! temporary logical unit
	196	INTEGER :: ji, jj, jk ! dummy loop indices
[479]	197	INTEGER :: ijloc
[508]	198	REAL(wp) :: zh1, zh2, zslope, zcst, zfcor ! temporary scalars
[2104]	199	REAL(wp) :: zt1 = 15._wp ! surface temperature value (EEL R5)
	200	REAL(wp) :: zt2 = 5._wp ! bottom temperature value (EEL R5)
	201	REAL(wp) :: zsal = 35.0_wp ! constant salinity (EEL R2, R5 and R6)
	202	REAL(wp) :: zueel = 0.1_wp ! constant uniform zonal velocity (EEL R5)
[508]	203	REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain
[3]	204	!!----------------------------------------------------------------------
	205
	206	SELECT CASE ( jp_cfg )
	207	! ! ====================
	208	CASE ( 5 ) ! EEL R5 configuration
	209	! ! ====================
[2104]	210	!
[3]	211	! set temperature field with a linear profile
	212	! -------------------------------------------
	213	IF(lwp) WRITE(numout,*)
	214	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
	215	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[2104]	216	!
[467]	217	zh1 = gdept_0( 1 )
	218	zh2 = gdept_0(jpkm1)
[2104]	219	!
[3]	220	zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
	221	zcst = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
[2104]	222	!
[3]	223	DO jk = 1, jpk
[2977]	224	tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
	225	tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
[3]	226	END DO
[2104]	227	!
[2977]	228	IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &
	229	& 1 , jpi , 5 , 1 , jpk , &
	230	& 1 , 1. , numout )
[2104]	231	!
[3]	232	! set salinity field to a constant value
	233	! --------------------------------------
	234	IF(lwp) WRITE(numout,*)
	235	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
	236	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[2104]	237	!
[2977]	238	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
	239	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
[2104]	240	!
[3]	241	! set the dynamics: U,V, hdiv, rot (and ssh if necessary)
	242	! ----------------
	243	! Start EEL5 configuration with barotropic geostrophic velocities
	244	! according the sshb and sshn SSH imposed.
[479]	245	! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
	246	! we use the Coriolis frequency at mid-channel.
	247	ub(:,:,:) = zueel * umask(:,:,:)
[3]	248	un(:,:,:) = ub(:,:,:)
[479]	249	ijloc = mj0(INT(jpjglo-1)/2)
	250	zfcor = ff(1,ijloc)
[2104]	251	!
[3]	252	DO jj = 1, jpjglo
[479]	253	zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)zueele1t(1,1)*zfcor/grav
[3]	254	END DO
[2104]	255	!
[479]	256	IF(lwp) THEN
	257	WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
	258	WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
	259	WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
	260	ENDIF
[2104]	261	!
[3]	262	DO jj = 1, nlcj
	263	DO ji = 1, nlci
	264	sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
	265	END DO
	266	END DO
	267	sshb(nlci+1:jpi, : ) = 0.e0 ! set to zero extra mpp columns
	268	sshb( : ,nlcj+1:jpj) = 0.e0 ! set to zero extra mpp rows
[2104]	269	!
[3]	270	sshn(:,:) = sshb(:,:) ! set now ssh to the before value
[2104]	271	!
[593]	272	IF( nn_rstssh /= 0 ) THEN
[2104]	273	nn_rstssh = 0 ! hand-made initilization of ssh
[593]	274	CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
[558]	275	ENDIF
[2104]	276	!
	277	CALL div_cur( nit000 ) ! horizontal divergence and relative vorticity (curl)
[3]	278	! N.B. the vertical velocity will be computed from the horizontal divergence field
	279	! in istate by a call to wzv routine
	280
	281
	282	! ! ==========================
	283	CASE ( 2 , 6 ) ! EEL R2 or R6 configuration
	284	! ! ==========================
[2104]	285	!
[3]	286	! set temperature field with a NetCDF file
	287	! ----------------------------------------
	288	IF(lwp) WRITE(numout,*)
	289	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
	290	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[2104]	291	!
[473]	292	CALL iom_open ( 'eel.initemp', inum )
[2977]	293	CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
[473]	294	CALL iom_close( inum )
[2104]	295	!
[2977]	296	tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb
[2104]	297	!
[2977]	298	IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &
	299	& 1 , jpi , 5 , 1 , jpk , &
	300	& 1 , 1. , numout )
[2104]	301	!
[3]	302	! set salinity field to a constant value
	303	! --------------------------------------
	304	IF(lwp) WRITE(numout,*)
	305	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
	306	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
[2104]	307	!
[2977]	308	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
	309	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
[2104]	310	!
[3]	311	! ! ===========================
	312	CASE DEFAULT ! NONE existing configuration
	313	! ! ===========================
[473]	314	WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
	315	CALL ctl_stop( ctmp1 )
[2104]	316	!
[3]	317	END SELECT
[2104]	318	!
[3]	319	END SUBROUTINE istate_eel
	320
	321
[93]	322	SUBROUTINE istate_gyre
	323	!!----------------------------------------------------------------------
	324	!! * ROUTINE istate_gyre *
	325	!!
	326	!! ** Purpose : Initialization of the dynamics and tracers for GYRE
	327	!! configuration (double gyre with rotated domain)
	328	!!
	329	!! ** Method : - set temprature field
	330	!! - set salinity field
	331	!!----------------------------------------------------------------------
[473]	332	INTEGER :: ji, jj, jk ! dummy loop indices
[508]	333	INTEGER :: inum ! temporary logical unit
	334	INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization
[93]	335	!!----------------------------------------------------------------------
	336
[434]	337	SELECT CASE ( ntsinit)
[93]	338
[434]	339	CASE ( 0 ) ! analytical T/S profil deduced from LEVITUS
	340	IF(lwp) WRITE(numout,*)
	341	IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
	342	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[93]	343
[434]	344	DO jk = 1, jpk
	345	DO jj = 1, jpj
	346	DO ji = 1, jpi
[2977]	347	tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 ) ) &
[434]	348	& * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2 &
	349	& + ( 15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) ) &
	350	& - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.) &
	351	& + 7. * (1500. - fsdept(ji,jj,jk)) / 1500. ) &
	352	& * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
[2977]	353	tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
	354	tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
[434]	355
[2977]	356	tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 ) ) &
[434]	357	& * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2 &
	358	& + ( 35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000. &
	359	& - 1.62 * TANH( (fsdept(ji,jj,jk) - 60. ) / 650. ) &
	360	& + 0.2 * TANH( (fsdept(ji,jj,jk) - 35. ) / 100. ) &
	361	& + 0.2 * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.) ) &
	362	& * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
[2977]	363	tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
	364	tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
[434]	365	END DO
[93]	366	END DO
	367	END DO
	368
[434]	369	CASE ( 1 ) ! T/S data fields read in dta_tem.nc/data_sal.nc files
	370	IF(lwp) WRITE(numout,*)
	371	IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
	372	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	373	IF(lwp) WRITE(numout,*) ' NetCDF FORMAT'
	374
	375	! Read temperature field
	376	! ----------------------
[473]	377	CALL iom_open ( 'data_tem', inum )
[2977]	378	CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) )
[473]	379	CALL iom_close( inum )
[434]	380
[2977]	381	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
	382	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
[434]	383
	384	! Read salinity field
	385	! -------------------
[473]	386	CALL iom_open ( 'data_sal', inum )
[2977]	387	CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) )
[473]	388	CALL iom_close( inum )
[434]	389
[2977]	390	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
	391	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
[434]	392
	393	END SELECT
	394
[93]	395	IF(lwp) THEN
	396	WRITE(numout,*)
	397	WRITE(numout,*) ' Initial temperature and salinity profiles:'
[467]	398	WRITE(numout, "(9x,' level gdept_0 temperature salinity ')" )
[2977]	399	WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_0(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )
[93]	400	ENDIF
	401
	402	END SUBROUTINE istate_gyre
	403
	404
[3]	405	SUBROUTINE istate_uvg
	406	!!----------------------------------------------------------------------
	407	!! * ROUTINE istate_uvg *
	408	!!
	409	!! ** Purpose : Compute the geostrophic velocities from (tn,sn) fields
	410	!!
	411	!! ** Method : Using the hydrostatic hypothesis the now hydrostatic
	412	!! pressure is computed by integrating the in-situ density from the
	413	!! surface to the bottom.
	414	!! p=integral [ rau*g dz ]
	415	!!----------------------------------------------------------------------
[359]	416	USE dynspg ! surface pressure gradient (dyn_spg routine)
[3]	417	USE divcur ! hor. divergence & rel. vorticity (div_cur routine)
	418	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	419
	420	INTEGER :: ji, jj, jk ! dummy loop indices
	421	INTEGER :: indic ! ???
[508]	422	REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars
[3162]	423	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
[3]	424	!!----------------------------------------------------------------------
[3162]	425	!
	426	CALL wrk_alloc( jpi, jpj, jpk, zprn)
	427	!
[3]	428	IF(lwp) WRITE(numout,*)
	429	IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
	430	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	431
	432	! Compute the now hydrostatic pressure
	433	! ------------------------------------
	434
[15]	435	zalfg = 0.5 * grav * rau0
[508]	436
	437	zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value
[3]	438
[508]	439	DO jk = 2, jpkm1 ! Vertical integration from the surface
[3]	440	zprn(:,:,jk) = zprn(:,:,jk-1) &
[359]	441	& + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
[3]	442	END DO
	443
	444	! Compute geostrophic balance
	445	! ---------------------------
	446	DO jk = 1, jpkm1
	447	DO jj = 2, jpjm1
	448	DO ji = fs_2, fs_jpim1 ! vertor opt.
	449	zmsv = 1. / MAX( umask(ji-1,jj+1,jk) + umask(ji ,jj+1,jk) &
	450	+ umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) , 1. )
	451	zphv = ( zprn(ji ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1) &
	452	+ ( zprn(ji+1,jj+1,jk) - zprn(ji ,jj+1,jk) ) * umask(ji ,jj+1,jk) / e1u(ji ,jj+1) &
	453	+ ( zprn(ji ,jj ,jk) - zprn(ji-1,jj ,jk) ) * umask(ji-1,jj ,jk) / e1u(ji-1,jj ) &
	454	+ ( zprn(ji+1,jj ,jk) - zprn(ji ,jj ,jk) ) * umask(ji ,jj ,jk) / e1u(ji ,jj )
	455	zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
	456
	457	zmsu = 1. / MAX( vmask(ji+1,jj ,jk) + vmask(ji ,jj ,jk) &
	458	+ vmask(ji+1,jj-1,jk) + vmask(ji ,jj-1,jk) , 1. )
	459	zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj ,jk) ) * vmask(ji+1,jj ,jk) / e2v(ji+1,jj ) &
	460	+ ( zprn(ji ,jj+1,jk) - zprn(ji ,jj ,jk) ) * vmask(ji ,jj ,jk) / e2v(ji ,jj ) &
	461	+ ( zprn(ji+1,jj ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1) &
	462	+ ( zprn(ji ,jj ,jk) - zprn(ji ,jj-1,jk) ) * vmask(ji ,jj-1,jk) / e2v(ji ,jj-1)
	463	zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
	464
	465	! Compute the geostrophic velocities
	466	un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji ,jj-1) )
	467	vn(ji,jj,jk) = 2. * zphv / ( ff(ji,jj) + ff(ji-1,jj ) )
	468	END DO
	469	END DO
	470	END DO
	471
	472	IF(lwp) WRITE(numout,*) ' we force to zero bottom velocity'
	473
	474	! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
	475	! to have a zero bottom velocity
	476
	477	DO jk = 1, jpkm1
	478	un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
	479	vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
	480	END DO
	481
	482	CALL lbc_lnk( un, 'U', -1. )
	483	CALL lbc_lnk( vn, 'V', -1. )
	484
	485	ub(:,:,:) = un(:,:,:)
	486	vb(:,:,:) = vn(:,:,:)
	487
	488	! WARNING !!!!!
	489	! after initializing u and v, we need to calculate the initial streamfunction bsf.
	490	! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
	491	! to do that, we call dyn_spg with a special trick:
[508]	492	! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
	493	! right value assuming the velocities have been set up in one time step.
	494	! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
	495	! sets up s false trend to calculate the barotropic streamfunction.
[3]	496
	497	ua(:,:,:) = ub(:,:,:) / rdt
	498	va(:,:,:) = vb(:,:,:) / rdt
	499
[359]	500	! calls dyn_spg. we assume euler time step, starting from rest.
[3]	501	indic = 0
[359]	502	CALL dyn_spg( nit000, indic ) ! surface pressure gradient
[3]	503
	504	! the new velocity is ua*rdt
	505
	506	CALL lbc_lnk( ua, 'U', -1. )
	507	CALL lbc_lnk( va, 'V', -1. )
	508
	509	ub(:,:,:) = ua(:,:,:) * rdt
	510	vb(:,:,:) = va(:,:,:) * rdt
	511	ua(:,:,:) = 0.e0
	512	va(:,:,:) = 0.e0
	513	un(:,:,:) = ub(:,:,:)
	514	vn(:,:,:) = vb(:,:,:)
	515
	516	! Compute the divergence and curl
	517
	518	CALL div_cur( nit000 ) ! now horizontal divergence and curl
	519
	520	hdivb(:,:,:) = hdivn(:,:,:) ! set the before to the now value
	521	rotb (:,:,:) = rotn (:,:,:) ! set the before to the now value
[508]	522	!
[3162]	523	CALL wrk_dealloc( jpi, jpj, jpk, zprn)
[2715]	524	!
[3]	525	END SUBROUTINE istate_uvg
	526
	527	!!=====================================================================
	528	END MODULE istate

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

Download in other formats: