Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

istate.F90 @ 11738

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 11738

Download in other formats: