Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

traadv_cen2.F90

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 19.1 KB

Rev	Line
[3]	1	MODULE traadv_cen2
[888]	2	!!======================================================================
	3	!! * MODULE traadv_cen2 *
[2528]	4	!! Ocean tracers: horizontal & vertical advective trend
[888]	5	!!======================================================================
[4990]	6	!! History : OPA ! 2001-08 (G. Madec, E. Durand) v8.2 trahad+trazad=traadv
	7	!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
	8	!! - ! 2004-08 (C. Talandier) New trends organization
[1559]	9	!! - ! 2005-11 (V. Garnier) Surface pressure gradient organization
	10	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) Step reorganization
	11	!! - ! 2006-07 (G. madec) add ups_orca_set routine
	12	!! 3.2 ! 2009-07 (G. Madec) add avmb, avtb in restart for cen2 advection
[2528]	13	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
[3]	14	!!----------------------------------------------------------------------
[503]	15
	16	!!----------------------------------------------------------------------
[2528]	17	!! tra_adv_cen2 : update the tracer trend with the advection trends using a 2nd order centered scheme
	18	!! ups_orca_set : allow mixed upstream/centered scheme in specific area (set for orca 2 and 4 only)
[3]	19	!!----------------------------------------------------------------------
[2528]	20	USE oce, ONLY: tsn ! now ocean temperature and salinity
[3]	21	USE dom_oce ! ocean space and time domain
[1037]	22	USE eosbn2 ! equation of state
[4990]	23	USE trd_oce ! trends: ocean variables
	24	USE trdtra ! trends manager: tracers
[888]	25	USE closea ! closed sea
	26	USE sbcrnf ! river runoffs
	27	USE in_out_manager ! I/O manager
[1537]	28	USE iom ! IOM library
[132]	29	USE diaptr ! poleward transport diagnostics
[1201]	30	USE zdf_oce ! ocean vertical physics
[2528]	31	USE trc_oce ! share passive tracers/Ocean variables
[2715]	32	USE lib_mpp ! MPP library
[3294]	33	USE wrk_nemo ! Memory Allocation
	34	USE timing ! Timing
[4990]	35	USE phycst
[3]	36
	37	IMPLICIT NONE
	38	PRIVATE
	39
[4990]	40	PUBLIC tra_adv_cen2 ! routine called by traadv.F90
[3]	41
[2715]	42	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
	43	! ! and in closed seas (orca 2 and 4 configurations)
[3]	44	!! * Substitutions
	45	# include "domzgr_substitute.h90"
	46	# include "vectopt_loop_substitute.h90"
	47	!!----------------------------------------------------------------------
[2528]	48	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	49	!! $Id$
[2528]	50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	51	!!----------------------------------------------------------------------
	52	CONTAINS
	53
[3294]	54	SUBROUTINE tra_adv_cen2( kt, kit000, cdtype, pun, pvn, pwn, &
[4990]	55	& ptb, ptn, pta, kjpt )
[3]	56	!!----------------------------------------------------------------------
	57	!! * ROUTINE tra_adv_cen2 *
	58	!!
	59	!! ** Purpose : Compute the now trend due to the advection of tracers
	60	!! and add it to the general trend of passive tracer equations.
	61	!!
	62	!! ** Method : The advection is evaluated by a second order centered
	63	!! scheme using now fields (leap-frog scheme). In specific areas
	64	!! (vicinity of major river mouths, some straits, or where tn is
[457]	65	!! approaching the freezing point) it is mixed with an upstream
[3]	66	!! scheme for stability reasons.
[457]	67	!! Part 0 : compute the upstream / centered flag
	68	!! (3D array, zind, defined at T-point (0<zind<1))
	69	!! Part I : horizontal advection
	70	!! * centered flux:
[2528]	71	!! zcenu = e2u*e3u un mi(ptn)
	72	!! zcenv = e1v*e3v vn mj(ptn)
[457]	73	!! * upstream flux:
[2528]	74	!! zupsu = e2u*e3u un (ptb(i) or ptb(i-1) ) [un>0 or <0]
	75	!! zupsv = e1v*e3v vn (ptb(j) or ptb(j-1) ) [vn>0 or <0]
[457]	76	!! * mixed upstream / centered horizontal advection scheme
[3]	77	!! zcofi = max(zind(i+1), zind(i))
	78	!! zcofj = max(zind(j+1), zind(j))
	79	!! zwx = zcofi * zupsu + (1-zcofi) * zcenu
	80	!! zwy = zcofj * zupsv + (1-zcofj) * zcenv
[457]	81	!! * horizontal advective trend (divergence of the fluxes)
[2528]	82	!! ztra = 1/(e1te2te3t) { di-1[zwx] + dj-1[zwy] }
[457]	83	!! * Add this trend now to the general trend of tracer (ta,sa):
[2528]	84	!! pta = pta + ztra
[4990]	85	!! * trend diagnostic (l_trdtra=T or l_trctra=T): the trend is
[457]	86	!! saved for diagnostics. The trends saved is expressed as
[4990]	87	!! Uh.gradh(T), i.e. save trend = ztra + ptn divn
[3]	88	!!
	89	!! Part II : vertical advection
	90	!! For temperature (idem for salinity) the advective trend is com-
	91	!! puted as follows :
[2528]	92	!! ztra = 1/e3t dk+1[ zwz ]
[3]	93	!! where the vertical advective flux, zwz, is given by :
	94	!! zwz = zcofk * zupst + (1-zcofk) * zcent
[457]	95	!! with
[2528]	96	!! zupsv = upstream flux = wn * (ptb(k) or ptb(k-1) ) [wn>0 or <0]
[3]	97	!! zcenu = centered flux = wn * mk(tn)
[457]	98	!! The surface boundary condition is :
[1528]	99	!! variable volume (lk_vvl = T) : zero advective flux
[2528]	100	!! lin. free-surf (lk_vvl = F) : wn(:,:,1) * ptn(:,:,1)
[3]	101	!! Add this trend now to the general trend of tracer (ta,sa):
[2528]	102	!! pta = pta + ztra
[4990]	103	!! Trend diagnostic (l_trdtra=T or l_trctra=T): the trend is
[457]	104	!! saved for diagnostics. The trends saved is expressed as :
[2528]	105	!! save trend = w.gradz(T) = ztra - ptn divn.
[3]	106	!!
[2528]	107	!! ** Action : - update pta with the now advective tracer trends
	108	!! - save trends if needed
[503]	109	!!----------------------------------------------------------------------
[2528]	110	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	111	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	112	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	113	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	114	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	115	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	116	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2715]	117	!
[4990]	118	INTEGER :: ji, jj, jk, jn, ikt ! dummy loop indices
	119	INTEGER :: ierr ! local integer
[2715]	120	REAL(wp) :: zbtr, ztra ! local scalars
	121	REAL(wp) :: zfp_ui, zfp_vj, zfp_w, zcofi ! - -
	122	REAL(wp) :: zfm_ui, zfm_vj, zfm_w, zcofj, zcofk ! - -
	123	REAL(wp) :: zupsut, zcenut, zupst ! - -
	124	REAL(wp) :: zupsvt, zcenvt, zcent, zice ! - -
[4990]	125	REAL(wp), POINTER, DIMENSION(:,:) :: zfzp, zpres ! 2D workspace
	126	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy ! 3D -
	127	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind ! - -
[3]	128	!!----------------------------------------------------------------------
[3294]	129	!
	130	IF( nn_timing == 1 ) CALL timing_start('tra_adv_cen2')
	131	!
[4990]	132	CALL wrk_alloc( jpi, jpj, zpres, zfzp )
	133	CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwz, zind )
[3294]	134	!
[3]	135
[3294]	136	IF( kt == kit000 ) THEN
[3]	137	IF(lwp) WRITE(numout,*)
[2528]	138	IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme on ', cdtype
	139	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
[3]	140	IF(lwp) WRITE(numout,*)
[11101]	141	IF(lwp .AND. lflush) CALL flush(numout)
[888]	142	!
[4990]	143	IF( .NOT. ALLOCATED( upsmsk ) ) THEN
[3294]	144	ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
	145	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array')
	146	ENDIF
	147
[2715]	148	!
	149	upsmsk(:,:) = 0._wp ! not upstream by default
[916]	150	!
[888]	151	IF( cp_cfg == "orca" ) CALL ups_orca_set ! set mixed Upstream/centered scheme near some straits
	152	! ! and in closed seas (orca2 and orca4 only)
[1537]	153	IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN ! Increase the background in the surface layers
	154	avmb(1) = 10. * avmb(1) ; avtb(1) = 10. * avtb(1)
	155	avmb(2) = 10. * avmb(2) ; avtb(2) = 10. * avtb(2)
	156	avmb(3) = 5. * avmb(3) ; avtb(3) = 5. * avtb(3)
	157	avmb(4) = 2.5 * avmb(4) ; avtb(4) = 2.5 * avtb(4)
	158	ENDIF
[3]	159	ENDIF
[2528]	160	!
[3]	161	! Upstream / centered scheme indicator
	162	! ------------------------------------
[1037]	163	!!gm not strickly exact : the freezing point should be computed at each ocean levels...
	164	!!gm not a big deal since cen2 is no more used in global ice-ocean simulations
[4990]	165	!!ch changes for ice shelf to retain standard behaviour elsewhere, even if not optimal
	166	DO jj = 1, jpj
	167	DO ji = 1, jpi
	168	ikt = mikt(ji,jj)
	169	IF (ikt > 1 ) THEN
	170	zpres(ji,jj) = grav * rau0 * fsdept(ji,jj,ikt) * 1.e-04
	171	ELSE
	172	zpres(ji,jj) = 0.0
	173	ENDIF
	174	END DO
	175	END DO
[6498]	176	CALL eos_fzp( tsn(:,:,1,jp_sal), zfzp(:,:), zpres(:,:) )
[3]	177	DO jk = 1, jpk
	178	DO jj = 1, jpj
	179	DO ji = 1, jpi
[1037]	180	! ! below ice covered area (if tn < "freezing"+0.1 )
[4990]	181	IF( tsn(ji,jj,jk,jp_tem) <= zfzp(ji,jj) + 0.1 ) THEN ; zice = 1._wp
	182	ELSE ; zice = 0._wp
[1037]	183	ENDIF
[888]	184	zind(ji,jj,jk) = MAX ( &
	185	rnfmsk(ji,jj) * rnfmsk_z(jk), & ! near runoff mouths (& closed sea outflows)
[1037]	186	upsmsk(ji,jj) , & ! some of some straits
	187	zice & ! below ice covered area (if tn < "freezing"+0.1 )
	188	& ) * tmask(ji,jj,jk)
[3]	189	END DO
	190	END DO
	191	END DO
	192
[2528]	193	DO jn = 1, kjpt
[503]	194	!
[2528]	195	! I. Horizontal advection
	196	! ====================
	197	!
[503]	198	DO jk = 1, jpkm1
[2528]	199	! ! Second order centered tracer flux at u- and v-points
	200	DO jj = 1, jpjm1
	201	!
	202	DO ji = 1, fs_jpim1 ! vector opt.
	203	! upstream indicator
	204	zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
	205	zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
[1559]	206	!
[2528]	207	! upstream scheme
	208	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
	209	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
	210	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
	211	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
	212	zupsut = zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn)
	213	zupsvt = zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn)
	214	! centered scheme
	215	zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) )
	216	zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) )
	217	! mixed centered / upstream scheme
	218	zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
	219	zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
[503]	220	END DO
	221	END DO
	222	END DO
[2528]	223
	224	! II. Vertical advection
	225	! ==================
[503]	226	!
[2528]	227	! ! Vertical advective fluxes
	228	zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero
	229	! ! Surface value :
	230	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable
[4990]	231	ELSE
	232	DO jj = 1, jpj ! vector opt.
	233	DO ji = 1, jpi ! vector opt.
	234	ikt = mikt(ji,jj)
	235	zwz(ji,jj,ikt ) = pwn(ji,jj,ikt) * ptn(ji,jj,ikt,jn) ! linear free surface
	236	zwz(ji,jj,1:ikt-1) = 0.e0
	237	END DO
	238	END DO
[2528]	239	ENDIF
	240	!
	241	DO jk = 2, jpk ! Second order centered tracer flux at w-point
[503]	242	DO jj = 2, jpjm1
	243	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	244	! upstream indicator
	245	zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
	246	! mixed centered / upstream scheme
	247	zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
	248	zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
	249	zupst = zfp_w * ptb(ji,jj,jk,jn) + zfm_w * ptb(ji,jj,jk-1,jn)
	250	! centered scheme
	251	zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )
	252	! mixed centered / upstream scheme
	253	zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
[503]	254	END DO
	255	END DO
	256	END DO
[216]	257
[2528]	258	! II. Divergence of advective fluxes
	259	! ----------------------------------
[503]	260	DO jk = 1, jpkm1
	261	DO jj = 2, jpjm1
	262	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	263	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	264	! advective trends
	265	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	266	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	267	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
	268	! advective trends added to the general tracer trends
	269	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[503]	270	END DO
	271	END DO
	272	END DO
[216]	273
[4990]	274	! ! trend diagnostics
	275	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. &
	276	&( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN
	277	CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
	278	CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
	279	CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
[2528]	280	END IF
	281	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[7179]	282	IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:) )
[2528]	283	!
[4990]	284	END DO
[2528]	285
[1537]	286	! --------------------------- required in restart file to ensure restartability)
	287	! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
[2528]	288	IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
[9321]	289	IF(nn_timing == 2) CALL timing_start('iom_rstput')
[1537]	290	CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
	291	CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
[9321]	292	IF(nn_timing == 2) CALL timing_stop('iom_rstput')
[1537]	293	ENDIF
[503]	294	!
[4990]	295	CALL wrk_dealloc( jpi, jpj, zpres, zfzp )
	296	CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwz, zind )
[2715]	297	!
[3294]	298	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_cen2')
	299	!
[3]	300	END SUBROUTINE tra_adv_cen2
[888]	301
	302
	303	SUBROUTINE ups_orca_set
	304	!!----------------------------------------------------------------------
	305	!! * ROUTINE ups_orca_set *
	306	!!
	307	!! ** Purpose : add a portion of upstream scheme in area where the
	308	!! centered scheme generates too strong overshoot
	309	!!
	310	!! ** Method : orca (R4 and R2) confiiguration setting. Set upsmsk
	311	!! array to nozero value in some straith.
	312	!!
	313	!! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
	314	!!----------------------------------------------------------------------
	315	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
	316	!!----------------------------------------------------------------------
[3294]	317	!
	318	IF( nn_timing == 1 ) CALL timing_start('ups_orca_set')
	319	!
[888]	320	! mixed upstream/centered scheme near river mouths
	321	! ------------------------------------------------
	322	SELECT CASE ( jp_cfg )
	323	! ! =======================
	324	CASE ( 4 ) ! ORCA_R4 configuration
	325	! ! =======================
	326	! ! Gibraltar Strait
	327	ii0 = 70 ; ii1 = 71
	328	ij0 = 52 ; ij1 = 53 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
	329	!
	330	! ! =======================
	331	CASE ( 2 ) ! ORCA_R2 configuration
	332	! ! =======================
	333	! ! Gibraltar Strait
	334	ij0 = 102 ; ij1 = 102
	335	ii0 = 138 ; ii1 = 138 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
	336	ii0 = 139 ; ii1 = 139 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
	337	ii0 = 140 ; ii1 = 140 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
	338	ij0 = 101 ; ij1 = 102
	339	ii0 = 141 ; ii1 = 141 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
	340	! ! Bab el Mandeb Strait
	341	ij0 = 87 ; ij1 = 88
	342	ii0 = 164 ; ii1 = 164 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
	343	ij0 = 88 ; ij1 = 88
	344	ii0 = 163 ; ii1 = 163 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
	345	ii0 = 162 ; ii1 = 162 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
	346	ii0 = 160 ; ii1 = 161 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
	347	ij0 = 89 ; ij1 = 89
	348	ii0 = 158 ; ii1 = 160 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
	349	ij0 = 90 ; ij1 = 90
	350	ii0 = 160 ; ii1 = 160 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
	351	! ! Sound Strait
	352	ij0 = 116 ; ij1 = 116
	353	ii0 = 144 ; ii1 = 144 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
	354	ii0 = 145 ; ii1 = 147 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
	355	ii0 = 148 ; ii1 = 148 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
	356	!
	357	END SELECT
	358
	359	! mixed upstream/centered scheme over closed seas
	360	! -----------------------------------------------
	361	CALL clo_ups( upsmsk(:,:) )
	362	!
[3294]	363	IF( nn_timing == 1 ) CALL timing_stop('ups_orca_set')
	364	!
[888]	365	END SUBROUTINE ups_orca_set
[3]	366
	367	!!======================================================================
	368	END MODULE traadv_cen2

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_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90

Download in other formats: