Context Navigation

trcdmp.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

File size: 18.3 KB

Rev	Line
[941]	1	MODULE trcdmp
	2	!!======================================================================
	3	!! * MODULE trcdmp *
	4	!! Ocean physics: internal restoring trend on passive tracers
	5	!!======================================================================
[2528]	6	!! History : OPA ! 1991-03 (O. Marti, G. Madec) Original code
	7	!! ! 1996-01 (G. Madec) statement function for e3
	8	!! ! 1997-05 (H. Loukos) adapted for passive tracers
	9	!! NEMO 9.0 ! 2004-03 (C. Ethe) free form + modules
	10	!! 3.2 ! 2007-02 (C. Deltel) Diagnose ML trends for passive tracers
	11	!! 3.3 ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC
[941]	12	!!----------------------------------------------------------------------
[4148]	13	#if defined key_top
[941]	14	!!----------------------------------------------------------------------
	15	!! trc_dmp : update the tracer trend with the internal damping
	16	!! trc_dmp_init : initialization, namlist read, parameters control
	17	!!----------------------------------------------------------------------
	18	USE oce_trc ! ocean dynamics and tracers variables
[1175]	19	USE trc ! ocean passive tracers variables
[941]	20	USE trcdta
[2528]	21	USE tradmp
	22	USE trdtra
[4990]	23	USE trd_oce
[9169]	24	!
[5102]	25	USE iom
[9169]	26	USE prtctl_trc ! Print control for debbuging
[941]	27
	28	IMPLICIT NONE
	29	PRIVATE
	30
[5836]	31	PUBLIC trc_dmp
	32	PUBLIC trc_dmp_clo
	33	PUBLIC trc_dmp_alloc
	34	PUBLIC trc_dmp_ini
[941]	35
[9169]	36	INTEGER , PUBLIC :: nn_zdmp_tr !: = 0/1/2 flag for damping in the mixed layer
	37	CHARACTER(LEN=200) , PUBLIC :: cn_resto_tr !: File containing restoration coefficient
[5836]	38
[4148]	39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: restotr ! restoring coeff. on tracers (s-1)
[2715]	40
[9169]	41	INTEGER, PARAMETER :: npncts = 8 ! number of closed sea
	42	INTEGER, DIMENSION(npncts) :: nctsi1, nctsj1 ! south-west closed sea limits (i,j)
	43	INTEGER, DIMENSION(npncts) :: nctsi2, nctsj2 ! north-east closed sea limits (i,j)
[941]	44
	45	!! * Substitutions
[5836]	46	# include "vectopt_loop_substitute.h90"
[12340]	47	# include "do_loop_substitute.h90"
[941]	48	!!----------------------------------------------------------------------
[10067]	49	!! NEMO/TOP 4.0 , NEMO Consortium (2018)
[5215]	50	!! $Id$
[10068]	51	!! Software governed by the CeCILL license (see ./LICENSE)
[941]	52	!!----------------------------------------------------------------------
	53	CONTAINS
	54
[2715]	55	INTEGER FUNCTION trc_dmp_alloc()
	56	!!----------------------------------------------------------------------
	57	!! * ROUTINE trc_dmp_alloc *
	58	!!----------------------------------------------------------------------
	59	ALLOCATE( restotr(jpi,jpj,jpk) , STAT=trc_dmp_alloc )
	60	!
	61	IF( trc_dmp_alloc /= 0 ) CALL ctl_warn('trc_dmp_alloc: failed to allocate array')
	62	!
	63	END FUNCTION trc_dmp_alloc
	64
	65
[11949]	66	SUBROUTINE trc_dmp( kt, Kbb, Kmm, ptr, Krhs )
[941]	67	!!----------------------------------------------------------------------
	68	!! * ROUTINE trc_dmp *
	69	!!
	70	!! ** Purpose : Compute the passive tracer trend due to a newtonian damping
	71	!! of the tracer field towards given data field and add it to the
	72	!! general tracer trends.
	73	!!
	74	!! ** Method : Newtonian damping towards trdta computed
	75	!! and add to the general tracer trends:
[11949]	76	!! tr(Kmm) = tr(Krhs) + restotr * (trdta - tr(Kbb))
[941]	77	!! The trend is computed either throughout the water column
	78	!! (nlmdmptr=0) or in area of weak vertical mixing (nlmdmptr=1) or
	79	!! below the well mixed layer (nlmdmptr=2)
	80	!!
[11949]	81	!! ** Action : - update the tracer trends tr(:,:,:,:,Krhs) with the newtonian
[941]	82	!! damping trends.
[4990]	83	!! - save the trends ('key_trdmxl_trc')
[1175]	84	!!----------------------------------------------------------------------
[11949]	85	INTEGER, INTENT(in ) :: kt ! ocean time-step index
	86	INTEGER, INTENT(in ) :: Kbb, Kmm, Krhs ! time level indices
	87	REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation
[6140]	88	!
	89	INTEGER :: ji, jj, jk, jn, jl ! dummy loop indices
	90	CHARACTER (len=22) :: charout
[9125]	91	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrd
	92	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrcdta ! 3D workspace
[941]	93	!!----------------------------------------------------------------------
[3294]	94	!
[9124]	95	IF( ln_timing ) CALL timing_start('trc_dmp')
[3294]	96	!
[9125]	97	IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) ! temporary save of trends
[4148]	98	!
	99	IF( nb_trcdta > 0 ) THEN ! Initialisation of tracer from a file that may also be used for damping
	100	!
[9125]	101	ALLOCATE( ztrcdta(jpi,jpj,jpk) ) ! Memory allocation
[4148]	102	! ! ===========
	103	DO jn = 1, jptra ! tracer loop
	104	! ! ===========
[11949]	105	IF( l_trdtrc ) ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs) ! save trends
[4148]	106	!
	107	IF( ln_trc_ini(jn) ) THEN ! update passive tracers arrays with input data read from file
[6140]	108	!
[4148]	109	jl = n_trc_index(jn)
[11949]	110	CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta ) ! read tracer data at nit000
[6140]	111	!
[4148]	112	SELECT CASE ( nn_zdmp_tr )
	113	!
	114	CASE( 0 ) !== newtonian damping throughout the water column ==!
[12340]	115	DO_3D_00_00( 1, jpkm1 )
	116	ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - ptr(ji,jj,jk,jn,Kbb) )
	117	END_3D
[6140]	118	!
[4148]	119	CASE ( 1 ) !== no damping in the turbocline (avt > 5 cm2/s) ==!
[12340]	120	DO_3D_00_00( 1, jpkm1 )
	121	IF( avt(ji,jj,jk) <= avt_c ) THEN
	122	ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - ptr(ji,jj,jk,jn,Kbb) )
	123	ENDIF
	124	END_3D
[6140]	125	!
[4148]	126	CASE ( 2 ) !== no damping in the mixed layer ==!
[12340]	127	DO_3D_00_00( 1, jpkm1 )
	128	IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
	129	ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - ptr(ji,jj,jk,jn,Kbb) )
	130	END IF
	131	END_3D
[6140]	132	!
[4148]	133	END SELECT
	134	!
	135	ENDIF
	136	!
	137	IF( l_trdtrc ) THEN
[11949]	138	ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs) - ztrtrd(:,:,:)
	139	CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_dmp, ztrtrd )
[4148]	140	END IF
	141	! ! ===========
	142	END DO ! tracer loop
	143	! ! ===========
[9125]	144	DEALLOCATE( ztrcdta )
[4148]	145	ENDIF
	146	!
[9125]	147	IF( l_trdtrc ) DEALLOCATE( ztrtrd )
[2528]	148	! ! print mean trends (used for debugging)
[12236]	149	IF( sn_cfctl%l_prttrc ) THEN
[6140]	150	WRITE(charout, FMT="('dmp ')")
	151	CALL prt_ctl_trc_info(charout)
[11949]	152	CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
[941]	153	ENDIF
[2528]	154	!
[9124]	155	IF( ln_timing ) CALL timing_stop('trc_dmp')
[3294]	156	!
[941]	157	END SUBROUTINE trc_dmp
	158
[6140]	159
[5836]	160	SUBROUTINE trc_dmp_ini
	161	!!----------------------------------------------------------------------
	162	!! * ROUTINE trc_dmp_ini *
	163	!!
	164	!! ** Purpose : Initialization for the newtonian damping
	165	!!
	166	!! ** Method : read the nammbf namelist and check the parameters
	167	!! called by trc_dmp at the first timestep (nittrc000)
	168	!!----------------------------------------------------------------------
[6140]	169	INTEGER :: ios, imask ! local integers
	170	!!
[5836]	171	NAMELIST/namtrc_dmp/ nn_zdmp_tr , cn_resto_tr
	172	!!----------------------------------------------------------------------
[6140]	173	!
[5836]	174	READ ( numnat_ref, namtrc_dmp, IOSTAT = ios, ERR = 909)
[11536]	175	909 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_dmp in reference namelist' )
[5836]	176	READ ( numnat_cfg, namtrc_dmp, IOSTAT = ios, ERR = 910)
[11536]	177	910 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtrc_dmp in configuration namelist' )
[5836]	178	IF(lwm) WRITE ( numont, namtrc_dmp )
	179
	180	IF(lwp) THEN ! Namelist print
	181	WRITE(numout,*)
	182	WRITE(numout,*) 'trc_dmp : Passive tracers newtonian damping'
	183	WRITE(numout,*) '~~~~~~~'
	184	WRITE(numout,*) ' Namelist namtrc_dmp : set damping parameter'
[9169]	185	WRITE(numout,*) ' mixed layer damping option nn_zdmp_tr = ', nn_zdmp_tr, '(zoom: forced to 0)'
	186	WRITE(numout,*) ' Restoration coeff file cn_resto_tr = ', cn_resto_tr
[5836]	187	ENDIF
[6701]	188	! ! Allocate arrays
	189	IF( trc_dmp_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trc_dmp_ini: unable to allocate arrays' )
[5836]	190	!
	191	SELECT CASE ( nn_zdmp_tr )
[9169]	192	CASE ( 0 ) ; IF(lwp) WRITE(numout,*) ' ===>> tracer damping throughout the water column'
	193	CASE ( 1 ) ; IF(lwp) WRITE(numout,*) ' ===>> no tracer damping in the turbocline (avt > 5 cm2/s)'
	194	CASE ( 2 ) ; IF(lwp) WRITE(numout,*) ' ===>> no tracer damping in the mixed layer'
[5836]	195	CASE DEFAULT
	196	WRITE(ctmp1,*) 'bad flag value for nn_zdmp_tr = ', nn_zdmp_tr
	197	CALL ctl_stop(ctmp1)
	198	END SELECT
	199
	200	IF( .NOT.lk_c1d ) THEN
[9169]	201	IF( .NOT.ln_tradmp ) &
	202	& CALL ctl_stop( 'passive tracer damping need ln_tradmp to compute damping coef.' )
[5836]	203	!
	204	! ! Read damping coefficients from file
	205	!Read in mask from file
	206	CALL iom_open ( cn_resto_tr, imask)
	207	CALL iom_get ( imask, jpdom_autoglo, 'resto', restotr)
	208	CALL iom_close( imask )
	209	!
	210	ENDIF
	211	!
	212	END SUBROUTINE trc_dmp_ini
	213
[6140]	214
[11949]	215	SUBROUTINE trc_dmp_clo( kt, Kbb, Kmm )
[4148]	216	!!---------------------------------------------------------------------
	217	!! * ROUTINE trc_dmp_clo *
	218	!!
	219	!! ** Purpose : Closed sea domain initialization
	220	!!
	221	!! ** Method : if a closed sea is located only in a model grid point
	222	!! we restore to initial data
	223	!!
	224	!! ** Action : nctsi1(), nctsj1() : south-west closed sea limits (i,j)
	225	!! nctsi2(), nctsj2() : north-east Closed sea limits (i,j)
	226	!!----------------------------------------------------------------------
[11949]	227	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	228	INTEGER, INTENT( in ) :: Kbb, Kmm ! time level indices
[4148]	229	!
[6701]	230	INTEGER :: ji , jj, jk, jn, jl, jc ! dummy loop indicesa
[6607]	231	INTEGER :: isrow ! local index
[6701]	232	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrcdta ! 3D workspace
[4148]	233	!!----------------------------------------------------------------------
[6607]	234
[4148]	235	IF( kt == nit000 ) THEN
	236	! initial values
	237	nctsi1(:) = 1 ; nctsi2(:) = 1
	238	nctsj1(:) = 1 ; nctsj2(:) = 1
	239
	240	! set the closed seas (in data domain indices)
	241	! -------------------
	242
[10213]	243	IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA") THEN
[4148]	244	!
[7646]	245	SELECT CASE ( nn_cfg )
[4148]	246	! ! =======================
[5506]	247	CASE ( 1 ) ! eORCA_R1 configuration
	248	! ! =======================
	249	isrow = 332 - jpjglo
	250	!
[9169]	251	nctsi1(1) = 333 ; nctsj1(1) = 243 - isrow ! Caspian Sea
[6607]	252	nctsi2(1) = 342 ; nctsj2(1) = 274 - isrow
[9169]	253	!
	254	nctsi1(2) = 198 ; nctsj1(2) = 258 - isrow ! Lake Superior
[6607]	255	nctsi2(2) = 204 ; nctsj2(2) = 262 - isrow
[9169]	256	!
	257	nctsi1(3) = 201 ; nctsj1(3) = 250 - isrow ! Lake Michigan
[6607]	258	nctsi2(3) = 203 ; nctsj2(3) = 256 - isrow
[9169]	259	!
	260	nctsi1(4) = 204 ; nctsj1(4) = 252 - isrow ! Lake Huron
[6607]	261	nctsi2(4) = 209 ; nctsj2(4) = 256 - isrow
[9169]	262	!
	263	nctsi1(5) = 206 ; nctsj1(5) = 249 - isrow ! Lake Erie
[6607]	264	nctsi2(5) = 209 ; nctsj2(5) = 251 - isrow
[9169]	265	!
	266	nctsi1(6) = 210 ; nctsj1(6) = 252 - isrow ! Lake Ontario
[6607]	267	nctsi2(6) = 212 ; nctsj2(6) = 252 - isrow
[9169]	268	!
	269	nctsi1(7) = 321 ; nctsj1(7) = 180 - isrow ! Victoria Lake
[6607]	270	nctsi2(7) = 322 ; nctsj2(7) = 189 - isrow
[9169]	271	!
	272	nctsi1(8) = 297 ; nctsj1(8) = 270 - isrow ! Baltic Sea
[6607]	273	nctsi2(8) = 308 ; nctsj2(8) = 293 - isrow
[6701]	274	!
[5506]	275	! ! =======================
[4148]	276	CASE ( 2 ) ! ORCA_R2 configuration
	277	! ! =======================
[9169]	278	!
	279	nctsi1(1) = 11 ; nctsj1(1) = 103 ! Caspian Sea
[4148]	280	nctsi2(1) = 17 ; nctsj2(1) = 112
[9169]	281	!
	282	nctsi1(2) = 97 ; nctsj1(2) = 107 ! Great North American Lakes
[4148]	283	nctsi2(2) = 103 ; nctsj2(2) = 111
[9169]	284	!
	285	nctsi1(3) = 174 ; nctsj1(3) = 107 ! Black Sea 1 : west part of the Black Sea
[4148]	286	nctsi2(3) = 181 ; nctsj2(3) = 112
[9169]	287	!
	288	nctsi1(4) = 2 ; nctsj1(4) = 107 ! Black Sea 2 : est part of the Black Sea
[4148]	289	nctsi2(4) = 6 ; nctsj2(4) = 112
[9169]	290	!
	291	nctsi1(5) = 145 ; nctsj1(5) = 116 ! Baltic Sea
[4148]	292	nctsi2(5) = 150 ; nctsj2(5) = 126
	293	! ! =======================
	294	CASE ( 4 ) ! ORCA_R4 configuration
	295	! ! =======================
[9169]	296	!
	297	nctsi1(1) = 4 ; nctsj1(1) = 53 ! Caspian Sea
[4148]	298	nctsi2(1) = 4 ; nctsj2(1) = 56
[9169]	299	!
	300	nctsi1(2) = 49 ; nctsj1(2) = 55 ! Great North American Lakes
[4148]	301	nctsi2(2) = 51 ; nctsj2(2) = 56
[9169]	302	!
	303	nctsi1(3) = 88 ; nctsj1(3) = 55 ! Black Sea
[4148]	304	nctsi2(3) = 91 ; nctsj2(3) = 56
[9169]	305	!
	306	nctsi1(4) = 75 ; nctsj1(4) = 59 ! Baltic Sea
[4148]	307	nctsi2(4) = 76 ; nctsj2(4) = 61
	308	! ! =======================
	309	CASE ( 025 ) ! ORCA_R025 configuration
	310	! ! =======================
[9169]	311	!
	312	nctsi1(1) = 1330 ; nctsj1(1) = 645 ! Caspian + Aral sea
[4148]	313	nctsi2(1) = 1400 ; nctsj2(1) = 795
[9169]	314	!
	315	nctsi1(2) = 1284 ; nctsj1(2) = 722 ! Azov Sea
[4148]	316	nctsi2(2) = 1304 ; nctsj2(2) = 747
	317	!
	318	END SELECT
	319	!
	320	ENDIF
	321	!
	322	! convert the position in local domain indices
	323	! --------------------------------------------
	324	DO jc = 1, npncts
	325	nctsi1(jc) = mi0( nctsi1(jc) )
	326	nctsj1(jc) = mj0( nctsj1(jc) )
[9169]	327	!
[4148]	328	nctsi2(jc) = mi1( nctsi2(jc) )
	329	nctsj2(jc) = mj1( nctsj2(jc) )
	330	END DO
	331	!
	332	ENDIF
	333
	334	! Restore close seas values to initial data
	335	IF( ln_trcdta .AND. nb_trcdta > 0 ) THEN ! Initialisation of tracer from a file that may also be used for damping
	336	!
	337	IF(lwp) WRITE(numout,*)
	338	IF(lwp) WRITE(numout,*) ' trc_dmp_clo : Restoring of nutrients on close seas at time-step kt = ', kt
	339	IF(lwp) WRITE(numout,*)
	340	!
[9125]	341	ALLOCATE( ztrcdta(jpi,jpj,jpk) ) ! Memory allocation
[6607]	342	!
[4148]	343	DO jn = 1, jptra
	344	IF( ln_trc_ini(jn) ) THEN ! update passive tracers arrays with input data read from file
	345	jl = n_trc_index(jn)
[11949]	346	CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta ) ! read tracer data at nit000
[4148]	347	DO jc = 1, npncts
	348	DO jk = 1, jpkm1
	349	DO jj = nctsj1(jc), nctsj2(jc)
	350	DO ji = nctsi1(jc), nctsi2(jc)
[11949]	351	tr(ji,jj,jk,jn,Kmm) = ztrcdta(ji,jj,jk)
	352	tr(ji,jj,jk,jn,Kbb) = tr(ji,jj,jk,jn,Kmm)
[9169]	353	END DO
	354	END DO
	355	END DO
	356	END DO
[4148]	357	ENDIF
[9169]	358	END DO
[9125]	359	DEALLOCATE( ztrcdta )
[4148]	360	ENDIF
	361	!
	362	END SUBROUTINE trc_dmp_clo
[6607]	363
[941]	364	#else
	365	!!----------------------------------------------------------------------
[4148]	366	!! Dummy module : No passive tracer
[941]	367	!!----------------------------------------------------------------------
	368	CONTAINS
	369	SUBROUTINE trc_dmp( kt ) ! Empty routine
	370	INTEGER, INTENT(in) :: kt
	371	WRITE(,) 'trc_dmp: You should not have seen this print! error?', kt
	372	END SUBROUTINE trc_dmp
	373	#endif
[4148]	374
[941]	375	!!======================================================================
	376	END MODULE trcdmp

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

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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/TRP/trcdmp.F90 @ 12340

Download in other formats: