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.

trcrad.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcrad.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
Property svn:keywords set to `Id`
File size: 14.7 KB

Rev	Line
[941]	1	MODULE trcrad
	2	!!======================================================================
	3	!! * MODULE trcrad *
	4	!! Ocean passive tracers: correction of negative concentrations
	5	!!======================================================================
	6	!! History : - ! 01-01 (O. Aumont & E. Kestenare) Original code
	7	!! 1.0 ! 04-03 (C. Ethe) free form F90
	8	!!----------------------------------------------------------------------
	9	#if defined key_top
	10	!!----------------------------------------------------------------------
	11	!! 'key_top' TOP models
	12	!!----------------------------------------------------------------------
	13	!! trc_rad : correction of negative concentrations
	14	!!----------------------------------------------------------------------
[9788]	15	USE par_trc ! need jptra, number of passive tracers
[941]	16	USE oce_trc ! ocean dynamics and tracers variables
[2528]	17	USE trc ! ocean passive tracers variables
[4990]	18	USE trd_oce
[2528]	19	USE trdtra
[941]	20	USE prtctl_trc ! Print control for debbuging
[10425]	21	USE lib_fortran
[941]	22
	23	IMPLICIT NONE
	24	PRIVATE
	25
[5836]	26	PUBLIC trc_rad
	27	PUBLIC trc_rad_ini
[941]	28
[5836]	29	LOGICAL , PUBLIC :: ln_trcrad !: flag to artificially correct negative concentrations
[10425]	30	REAL(wp), DIMENSION(:,:), ALLOCATABLE:: gainmass
[5836]	31
[941]	32	!!----------------------------------------------------------------------
[10067]	33	!! NEMO/TOP 4.0 , NEMO Consortium (2018)
[7753]	34	!! $Id$
[10068]	35	!! Software governed by the CeCILL license (see ./LICENSE)
[941]	36	!!----------------------------------------------------------------------
	37	CONTAINS
	38
[10946]	39	SUBROUTINE trc_rad( kt, Kmm, Krhs )
[941]	40	!!----------------------------------------------------------------------
	41	!! * ROUTINE trc_rad *
	42	!!
	43	!! ** Purpose : "crappy" routine to correct artificial negative
	44	!! concentrations due to isopycnal scheme
	45	!!
	46	!! ** Method : - PISCES or LOBSTER: Set negative concentrations to zero
	47	!! while computing the corresponding tracer content that
	48	!! is added to the tracers. Then, adjust the tracer
	49	!! concentration using a multiplicative factor so that
	50	!! the total tracer concentration is preserved.
	51	!! - CFC: simply set to zero the negative CFC concentration
	52	!! (the total CFC content is not strictly preserved)
	53	!!----------------------------------------------------------------------
[9169]	54	INTEGER, INTENT(in) :: kt ! ocean time-step index
[10946]	55	INTEGER, INTENT(in) :: Kmm, Krhs ! time level indices
[9169]	56	!
[941]	57	CHARACTER (len=22) :: charout
	58	!!----------------------------------------------------------------------
[3294]	59	!
[9124]	60	IF( ln_timing ) CALL timing_start('trc_rad')
[3294]	61	!
[10946]	62	IF( ln_age ) CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_age , jp_age ) ! AGE
	63	IF( ll_cfc ) CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_cfc0, jp_cfc1 ) ! CFC model
	64	IF( ln_c14 ) CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_c14 , jp_c14 ) ! C14
	65	IF( ln_pisces ) CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_pcs0, jp_pcs1, cpreserv='Y' ) ! PISCES model
	66	IF( ln_my_trc ) CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_myt0, jp_myt1 ) ! MY_TRC model
[1003]	67	!
	68	IF(ln_ctl) THEN ! print mean trends (used for debugging)
	69	WRITE(charout, FMT="('rad')")
	70	CALL prt_ctl_trc_info( charout )
	71	CALL prt_ctl_trc( tab4d=trn, mask=tmask, clinfo=ctrcnm )
	72	ENDIF
	73	!
[9124]	74	IF( ln_timing ) CALL timing_stop('trc_rad')
[3294]	75	!
[1003]	76	END SUBROUTINE trc_rad
	77
[9169]	78
[5836]	79	SUBROUTINE trc_rad_ini
	80	!!---------------------------------------------------------------------
	81	!! * ROUTINE trc _rad_ini *
	82	!!
[9169]	83	!! ** Purpose : read namelist options
[5836]	84	!!----------------------------------------------------------------------
[9169]	85	INTEGER :: ios ! Local integer output status for namelist read
	86	!!
[5836]	87	NAMELIST/namtrc_rad/ ln_trcrad
	88	!!----------------------------------------------------------------------
	89	!
	90	REWIND( numnat_ref ) ! namtrc_rad in reference namelist
	91	READ ( numnat_ref, namtrc_rad, IOSTAT = ios, ERR = 907)
[9169]	92	907 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_rad in reference namelist', lwp )
[5836]	93	REWIND( numnat_cfg ) ! namtrc_rad in configuration namelist
	94	READ ( numnat_cfg, namtrc_rad, IOSTAT = ios, ERR = 908 )
[9169]	95	908 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtrc_rad in configuration namelist', lwp )
	96	IF(lwm) WRITE( numont, namtrc_rad )
[5836]	97
	98	IF(lwp) THEN ! ! Control print
	99	WRITE(numout,*)
[9169]	100	WRITE(numout,*) 'trc_rad : Correct artificial negative concentrations '
	101	WRITE(numout,*) '~~~~~~~ '
[5836]	102	WRITE(numout,*) ' Namelist namtrc_rad : treatment of negative concentrations'
[9169]	103	WRITE(numout,*) ' correct artificially negative concen. or not ln_trcrad = ', ln_trcrad
	104	WRITE(numout,*)
	105	IF( ln_trcrad ) THEN ; WRITE(numout,*) ' ===>> ensure the global tracer conservation'
	106	ELSE ; WRITE(numout,*) ' ===>> NO strict global tracer conservation'
	107	ENDIF
[5836]	108	ENDIF
	109	!
[10425]	110	ALLOCATE( gainmass(jptra,2) )
	111	gainmass(:,:) = 0.
	112	!
[5836]	113	END SUBROUTINE trc_rad_ini
	114
[9169]	115
[10946]	116	SUBROUTINE trc_rad_sms( kt, Kmm, Krhs, ptrb, ptrn, jp_sms0, jp_sms1, cpreserv )
[1003]	117	!!-----------------------------------------------------------------------------
	118	!! * ROUTINE trc_rad_sms *
	119	!!
	120	!! ** Purpose : "crappy" routine to correct artificial negative
	121	!! concentrations due to isopycnal scheme
	122	!!
	123	!! ** Method : 2 cases :
	124	!! - Set negative concentrations to zero while computing
	125	!! the corresponding tracer content that is added to the
	126	!! tracers. Then, adjust the tracer concentration using
	127	!! a multiplicative factor so that the total tracer
	128	!! concentration is preserved.
	129	!! - simply set to zero the negative CFC concentration
	130	!! (the total content of concentration is not strictly preserved)
	131	!!--------------------------------------------------------------------------------
[9169]	132	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[10946]	133	INTEGER , INTENT(in ) :: Kmm, Krhs ! time level indices
[9169]	134	INTEGER , INTENT(in ) :: jp_sms0, jp_sms1 ! First & last index of the passive tracer model
	135	REAL(wp), DIMENSION (jpi,jpj,jpk,jptra), INTENT(inout) :: ptrb , ptrn ! before and now traceur concentration
	136	CHARACTER( len = 1), OPTIONAL , INTENT(in ) :: cpreserv ! flag to preserve content or not
	137	!
[10425]	138	INTEGER :: ji, ji2, jj, jj2, jk, jn ! dummy loop indices
	139	INTEGER :: icnt
	140	LOGICAL :: lldebug = .FALSE. ! local logical
	141	REAL(wp):: zcoef, zs2rdt, ztotmass
	142	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrneg, ztrpos
	143	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrd ! workspace arrays
[1003]	144	!!----------------------------------------------------------------------
[9169]	145	!
[10425]	146	IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) )
	147	zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) )
[9169]	148	!
	149	IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==!
	150	!
[10425]	151	ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) )
	152
[1003]	153	DO jn = jp_sms0, jp_sms1
[10425]	154	ztrneg(:,:,jn) = SUM( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values
	155	ztrpos(:,:,jn) = SUM( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values
	156	END DO
	157	CALL sum3x3( ztrneg )
	158	CALL sum3x3( ztrpos )
	159
	160	DO jn = jp_sms0, jp_sms1
[9169]	161	!
[10425]	162	IF( l_trdtrc ) ztrtrd(:,:,:) = ptrb(:,:,:,jn) ! save input trb for trend computation
[9169]	163	!
[10425]	164	DO jk = 1, jpkm1
	165	DO jj = 1, jpj
	166	DO ji = 1, jpi
	167	IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box
	168	!
	169	ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * tmask(ji,jj,jk) ! really needed?
	170	IF( ptrb(ji,jj,jk,jn) < 0. ) ptrb(ji,jj,jk,jn) = 0. ! supress negative values
	171	IF( ptrb(ji,jj,jk,jn) > 0. ) THEN ! use positive values to compensate mass gain
	172	zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptrb > 0
	173	ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef
	174	IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value
	175	gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk) ! we are adding mass...
	176	ptrb(ji,jj,jk,jn) = 0. ! limit the compensation to keep positive value
	177	ENDIF
	178	ENDIF
	179	!
	180	ENDIF
	181	END DO
	182	END DO
	183	END DO
	184	!
[3294]	185	IF( l_trdtrc ) THEN
[10425]	186	ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
[10946]	187	CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling
[3294]	188	ENDIF
[9169]	189	!
[10425]	190	END DO
	191
	192	IF( kt == nitend ) THEN
	193	CALL mpp_sum( 'trcrad', gainmass(:,1) )
	194	DO jn = jp_sms0, jp_sms1
	195	IF( gainmass(jn,1) > 0. ) THEN
	196	ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) )
	197	IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn &
	198	& , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1)
	199	END IF
	200	END DO
	201	ENDIF
	202
	203	DO jn = jp_sms0, jp_sms1
	204	ztrneg(:,:,jn) = SUM( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values
	205	ztrpos(:,:,jn) = SUM( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values
	206	END DO
	207	CALL sum3x3( ztrneg )
	208	CALL sum3x3( ztrpos )
	209
	210	DO jn = jp_sms0, jp_sms1
[9169]	211	!
[10425]	212	IF( l_trdtrc ) ztrtrd(:,:,:) = ptrn(:,:,:,jn) ! save input trb for trend computation
[9169]	213	!
[10425]	214	DO jk = 1, jpkm1
	215	DO jj = 1, jpj
	216	DO ji = 1, jpi
	217	IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box
	218	!
	219	ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * tmask(ji,jj,jk) ! really needed?
	220	IF( ptrn(ji,jj,jk,jn) < 0. ) ptrn(ji,jj,jk,jn) = 0. ! supress negative values
	221	IF( ptrn(ji,jj,jk,jn) > 0. ) THEN ! use positive values to compensate mass gain
	222	zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptrb > 0
	223	ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef
	224	IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value
	225	gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk) ! we are adding mass...
	226	ptrn(ji,jj,jk,jn) = 0. ! limit the compensation to keep positive value
	227	ENDIF
	228	ENDIF
	229	!
	230	ENDIF
	231	END DO
[941]	232	END DO
[10425]	233	END DO
[941]	234	!
[1175]	235	IF( l_trdtrc ) THEN
[10425]	236	ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
[10946]	237	CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd ) ! standard trend handling
[1175]	238	ENDIF
[9169]	239	!
[941]	240	END DO
[10425]	241
	242	IF( kt == nitend ) THEN
	243	CALL mpp_sum( 'trcrad', gainmass(:,2) )
	244	DO jn = jp_sms0, jp_sms1
	245	IF( gainmass(jn,2) > 0. ) THEN
	246	ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) )
	247	WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn &
	248	& , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1)
	249	END IF
	250	END DO
	251	ENDIF
	252
	253	DEALLOCATE( ztrneg, ztrpos )
[941]	254	!
[9169]	255	ELSE !== total CFC content is NOT strictly preserved ==!
[1003]	256	!
	257	DO jn = jp_sms0, jp_sms1
[9169]	258	!
[10425]	259	IF( l_trdtrc ) ztrtrd(:,:,:) = ptrb(:,:,:,jn) ! save input trb for trend computation
	260	!
	261	WHERE( ptrb(:,:,:,jn) < 0. ) ptrb(:,:,:,jn) = 0.
	262	!
[9169]	263	IF( l_trdtrc ) THEN
[10425]	264	ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
[10946]	265	CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling
[9169]	266	ENDIF
	267	!
[10425]	268	IF( l_trdtrc ) ztrtrd(:,:,:) = ptrn(:,:,:,jn) ! save input trn for trend computation
[9169]	269	!
[10425]	270	WHERE( ptrn(:,:,:,jn) < 0. ) ptrn(:,:,:,jn) = 0.
	271	!
[7753]	272	IF( l_trdtrc ) THEN
[10425]	273	ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
[10946]	274	CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd ) ! standard trend handling
[1257]	275	ENDIF
	276	!
[9169]	277	END DO
	278	!
[941]	279	ENDIF
[9169]	280	!
[10425]	281	IF( l_trdtrc ) DEALLOCATE( ztrtrd )
[9169]	282	!
	283	END SUBROUTINE trc_rad_sms
[1175]	284
[941]	285	#else
	286	!!----------------------------------------------------------------------
	287	!! Dummy module : NO TOP model
	288	!!----------------------------------------------------------------------
	289	CONTAINS
	290	SUBROUTINE trc_rad( kt ) ! Empty routine
	291	INTEGER, INTENT(in) :: kt
	292	WRITE(,) 'trc_rad: You should not have seen this print! error?', kt
	293	END SUBROUTINE trc_rad
	294	#endif
	295
	296	!!======================================================================
	297	END MODULE trcrad

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

Download in other formats: