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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
Property svn:keywords set to `Id`
File size: 10.0 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
[7646]	16	!! 3.7 ! 2016-04 (S. Flavoni) introduce user defined initial state
[508]	17	!!----------------------------------------------------------------------
[3]	18
	19	!!----------------------------------------------------------------------
	20	!! istate_init : initial state setting
	21	!! istate_uvg : initial velocity in geostropic balance
	22	!!----------------------------------------------------------------------
[7646]	23	USE oce ! ocean dynamics and active tracers
	24	USE dom_oce ! ocean space and time domain
	25	USE daymod ! calendar
	26	USE divhor ! horizontal divergence (div_hor routine)
	27	USE dtatsd ! data temperature and salinity (dta_tsd routine)
	28	USE dtauvd ! data: U & V current (dta_uvd routine)
[593]	29	USE domvvl ! varying vertical mesh
[6140]	30	USE iscplrst ! ice sheet coupling
[7646]	31	USE wet_dry ! wetting and drying (needed for wad_istate)
	32	USE usrdef_istate ! User defined initial state
[4990]	33	!
	34	USE in_out_manager ! I/O manager
	35	USE iom ! I/O library
[2715]	36	USE lib_mpp ! MPP library
[3680]	37	USE restart ! restart
[3294]	38	USE wrk_nemo ! Memory allocation
	39	USE timing ! Timing
[2715]	40
[3]	41	IMPLICIT NONE
	42	PRIVATE
	43
[508]	44	PUBLIC istate_init ! routine called by step.F90
[3]	45
	46	!! * Substitutions
	47	# include "vectopt_loop_substitute.h90"
	48	!!----------------------------------------------------------------------
[4990]	49	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
[888]	50	!! $Id$
[2715]	51	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	52	!!----------------------------------------------------------------------
	53	CONTAINS
	54
	55	SUBROUTINE istate_init
	56	!!----------------------------------------------------------------------
	57	!! * ROUTINE istate_init *
	58	!!
[508]	59	!! ** Purpose : Initialization of the dynamics and tracer fields.
[3]	60	!!----------------------------------------------------------------------
[7698]	61	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[5163]	62	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zuvd ! U & V data workspace
[3294]	63	!!----------------------------------------------------------------------
	64	!
[4990]	65	IF( nn_timing == 1 ) CALL timing_start('istate_init')
[3294]	66	!
[5836]	67	IF(lwp) WRITE(numout,*)
[7646]	68	IF(lwp) WRITE(numout,*) 'istate_init : Initialization of the dynamics and tracers'
	69	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[3]	70
[7646]	71	!!gm Why not include in the first call of dta_tsd ?
	72	!!gm probably associated with the use of internal damping...
[5836]	73	CALL dta_tsd_init ! Initialisation of T & S input data
[7646]	74	!!gm to be moved in usrdef of C1D case
	75	! IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data
	76	!!gm
[7698]	77	!$OMP PARALLEL
	78	DO jn = 1, jpts
	79	!$OMP DO schedule(static) private(jk, jj, ji)
	80	DO jk = 1, jpk
	81	DO jj = 1, jpj
	82	DO ji = 1, jpi
	83	tsa (ji,jj,jk,jn) = 0._wp ! set one for all to 0 at level jpk
	84	rab_b(ji,jj,jk,jn) = 0._wp ; rab_n(ji,jj,jk,jn) = 0._wp ! set one for all to 0 at level jpk
	85	END DO
	86	END DO
	87	END DO
	88	END DO
	89	!$OMP DO schedule(static) private(jk, jj, ji)
	90	DO jk = 1, jpk
	91	DO jj = 1, jpj
	92	DO ji = 1, jpi
	93	rhd (ji,jj,jk ) = 0._wp ; rhop (ji,jj,jk ) = 0._wp ! set one for all to 0 at level jpk
	94	rn2b (ji,jj,jk ) = 0._wp ; rn2 (ji,jj,jk ) = 0._wp ! set one for all to 0 at levels 1 and jpk
	95	END DO
	96	END DO
	97	END DO
	98	!$OMP END PARALLEL
[3]	99
[15]	100	IF( ln_rstart ) THEN ! Restart from a file
[3]	101	! ! -------------------
[6140]	102	CALL rst_read ! Read the restart file
[7646]	103	IF (ln_iscpl) CALL iscpl_stp ! extrapolate restart to wet and dry
[6140]	104	CALL day_init ! model calendar (using both namelist and restart infos)
[7646]	105	!
	106	ELSE ! Start from rest
[3]	107	! ! ---------------
[7646]	108	numror = 0 ! define numror = 0 -> no restart file to read
	109	neuler = 0 ! Set time-step indicator at nit000 (euler forward)
	110	CALL day_init ! model calendar (using both namelist and restart infos)
	111	! ! Initialization of ocean to zero
[508]	112	!
[7646]	113	IF( ln_tsd_init ) THEN
	114	CALL dta_tsd( nit000, tsb ) ! read 3D T and S data at nit000
	115	!
[7698]	116	!$OMP PARALLEL
	117	!$OMP DO schedule(static) private(jj, ji)
	118	DO jj = 1, jpj
	119	DO ji = 1, jpi
	120	sshb (ji,jj) = 0._wp ! set the ocean at rest
	121	END DO
	122	END DO
	123	!$OMP END DO NOWAIT
	124	!$OMP DO schedule(static) private(jk, jj, ji)
	125	DO jk = 1, jpk
	126	DO jj = 1, jpj
	127	DO ji = 1, jpi
	128	ub (ji,jj,jk) = 0._wp
	129	vb (ji,jj,jk) = 0._wp
	130	END DO
	131	END DO
	132	END DO
	133	!$OMP END PARALLEL
[7646]	134	!
	135	ELSE ! user defined initial T and S
	136	CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb )
[3]	137	ENDIF
[7698]	138	!$OMP PARALLEL
	139	DO jn = 1, jpts
	140	!$OMP DO schedule(static) private(jk, jj, ji)
	141	DO jk = 1, jpk
	142	DO jj = 1, jpj
	143	DO ji = 1, jpi
	144	tsn (ji,jj,jk,jn) = tsb (ji,jj,jk,jn) ! set now values from to before ones
	145	END DO
	146	END DO
	147	END DO
	148	END DO
	149	!$OMP DO schedule(static) private(jk, jj, ji)
	150	DO jk = 1, jpk
	151	DO jj = 1, jpj
	152	DO ji = 1, jpi
	153	un (ji,jj,jk) = ub (ji,jj,jk)
	154	vn (ji,jj,jk) = vb (ji,jj,jk)
	155	END DO
	156	END DO
	157	END DO
	158	!$OMP DO schedule(static) private(jj, ji)
	159	DO jj = 1, jpj
	160	DO ji = 1, jpi
	161	sshn (ji,jj) = sshb(ji,jj)
	162	hdivn(ji,jj,jpk) = 0._wp ! bottom divergence set one for 0 to zero at jpk level
	163	END DO
	164	END DO
	165	!$OMP END PARALLEL
[7646]	166	CALL div_hor( 0 ) ! compute interior hdivn value
	167	!!gm hdivn(:,:,:) = 0._wp
	168
	169	!!gm POTENTIAL BUG :
	170	!!gm ISSUE : if sshb /= 0 then, in non linear free surface, the e3._n, e3._b should be recomputed
	171	!! as well as gdept and gdepw.... !!!!!
	172	!! ===>>>> probably a call to domvvl initialisation here....
	173
	174
[6140]	175	!
[7646]	176	!!gm to be moved in usrdef of C1D case
	177	! IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
	178	! CALL wrk_alloc( jpi,jpj,jpk,2, zuvd )
	179	! CALL dta_uvd( nit000, zuvd )
	180	! ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:)
	181	! vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:)
	182	! CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd )
	183	! ENDIF
	184	!
[6140]	185	!!gm This is to be changed !!!!
[7646]	186	! ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
	187	! IF( .NOT.ln_linssh ) THEN
	188	! DO jk = 1, jpk
	189	! e3t_b(:,:,jk) = e3t_n(:,:,jk)
	190	! END DO
	191	! ENDIF
[6140]	192	!!gm
[2148]	193	!
[7646]	194	ENDIF
[4354]	195	!
	196	! Initialize "now" and "before" barotropic velocities:
[6140]	197	! Do it whatever the free surface method, these arrays being eventually used
[4354]	198	!
[7698]	199	!$OMP PARALLEL
	200	!$OMP DO schedule(static) private(jj, ji)
	201	DO jj = 1, jpj
	202	DO ji = 1, jpi
	203	un_b(ji,jj) = 0._wp ; vn_b(ji,jj) = 0._wp
	204	ub_b(ji,jj) = 0._wp ; vb_b(ji,jj) = 0._wp
	205	END DO
	206	END DO
[4354]	207	!
[7646]	208	!!gm the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
[4354]	209	DO jk = 1, jpkm1
[7698]	210	!$OMP DO schedule(static) private(jj, ji)
[4354]	211	DO jj = 1, jpj
	212	DO ji = 1, jpi
[6140]	213	un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
	214	vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
[4354]	215	!
[6140]	216	ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
	217	vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
[4354]	218	END DO
	219	END DO
	220	END DO
	221	!
[7698]	222	!$OMP DO schedule(static) private(jj, ji)
	223	DO jj = 1, jpj
	224	DO ji = 1, jpi
	225	un_b(ji,jj) = un_b(ji,jj) * r1_hu_n(ji,jj)
	226	vn_b(ji,jj) = vn_b(ji,jj) * r1_hv_n(ji,jj)
	227	!
	228	ub_b(ji,jj) = ub_b(ji,jj) * r1_hu_b(ji,jj)
	229	vb_b(ji,jj) = vb_b(ji,jj) * r1_hv_b(ji,jj)
	230	END DO
	231	END DO
	232	!$OMP END PARALLEL
[4354]	233	!
[4990]	234	IF( nn_timing == 1 ) CALL timing_stop('istate_init')
[3294]	235	!
[3]	236	END SUBROUTINE istate_init
	237
[7646]	238	!!======================================================================
[3]	239	END MODULE istate

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

Download in other formats: