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

Context Navigation

source: branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 2807

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

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

Download in other formats: