Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

diaar5.F90

Last change on this file was 11035, checked in by frrh, 5 years ago
Apply changes from UKMO GMED ticket 453 to fix OOB bugs in the calculation of diagnostic tnpeo and to apply the correct calculation as per the current NEMO vn4.0 and 3.6_STABLE code bases, in place of the existing incorrect formula. Thus, tnpeo diagnostic values will change following this update. Model prognostics and eveolution are unaffected.
File size: 13.9 KB

Rev	Line
[1756]	1	MODULE diaar5
	2	!!======================================================================
	3	!! * MODULE diaar5 *
	4	!! AR5 diagnostics
	5	!!======================================================================
[2528]	6	!! History : 3.2 ! 2009-11 (S. Masson) Original code
	7	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
[1756]	8	!!----------------------------------------------------------------------
[2715]	9	#if defined key_diaar5 \|\| defined key_esopa
[1756]	10	!!----------------------------------------------------------------------
	11	!! 'key_diaar5' : activate ar5 diagnotics
	12	!!----------------------------------------------------------------------
[2528]	13	!! dia_ar5 : AR5 diagnostics
	14	!! dia_ar5_init : initialisation of AR5 diagnostics
[1756]	15	!!----------------------------------------------------------------------
	16	USE oce ! ocean dynamics and active tracers
	17	USE dom_oce ! ocean space and time domain
[2528]	18	USE eosbn2 ! equation of state (eos_bn2 routine)
[1756]	19	USE lib_mpp ! distribued memory computing library
	20	USE iom ! I/O manager library
[3294]	21	USE timing ! preformance summary
	22	USE wrk_nemo ! working arrays
[5253]	23	USE fldread ! type FLD_N
	24	USE phycst ! physical constant
	25	USE in_out_manager ! I/O manager
[7179]	26	USE zdfddm
	27	USE zdf_oce
[1756]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[2528]	32	PUBLIC dia_ar5 ! routine called in step.F90 module
	33	PUBLIC dia_ar5_init ! routine called in opa.F90 module
[2715]	34	PUBLIC dia_ar5_alloc ! routine called in nemogcm.F90 module
[1756]	35
	36	LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .TRUE. ! coupled flag
	37
[2528]	38	REAL(wp) :: vol0 ! ocean volume (interior domain)
	39	REAL(wp) :: area_tot ! total ocean surface (interior domain)
[2715]	40	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: area ! cell surface (interior domain)
	41	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: thick0 ! ocean thickness (interior domain)
	42	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sn0 ! initial salinity
[1756]	43
	44	!! * Substitutions
	45	# include "domzgr_substitute.h90"
[7179]	46	# include "zdfddm_substitute.h90"
[1756]	47	!!----------------------------------------------------------------------
[2528]	48	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	49	!! $Id$
	50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1756]	51	!!----------------------------------------------------------------------
	52	CONTAINS
	53
[2715]	54	FUNCTION dia_ar5_alloc()
	55	!!----------------------------------------------------------------------
	56	!! * ROUTINE dia_ar5_alloc *
	57	!!----------------------------------------------------------------------
	58	INTEGER :: dia_ar5_alloc
	59	!!----------------------------------------------------------------------
	60	!
	61	ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )
	62	!
	63	IF( lk_mpp ) CALL mpp_sum ( dia_ar5_alloc )
	64	IF( dia_ar5_alloc /= 0 ) CALL ctl_warn('dia_ar5_alloc: failed to allocate arrays')
	65	!
	66	END FUNCTION dia_ar5_alloc
	67
	68
[1756]	69	SUBROUTINE dia_ar5( kt )
	70	!!----------------------------------------------------------------------
	71	!! * ROUTINE dia_ar5 *
	72	!!
[2528]	73	!! ** Purpose : compute and output some AR5 diagnostics
[1756]	74	!!----------------------------------------------------------------------
[2715]	75	!
[1756]	76	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	77	!
[1756]	78	INTEGER :: ji, jj, jk ! dummy loop arguments
	79	REAL(wp) :: zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass
[7179]	80	REAL(wp) :: zaw, zbw, zrw
[3294]	81	!
	82	REAL(wp), POINTER, DIMENSION(:,:) :: zarea_ssh , zbotpres ! 2D workspace
[11035]	83	REAL(wp), POINTER, DIMENSION(:,:) :: zpe ! 2D workspace
[3294]	84	REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhd , zrhop ! 3D workspace
	85	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztsn ! 4D workspace
[1756]	86	!!--------------------------------------------------------------------
[3294]	87	IF( nn_timing == 1 ) CALL timing_start('dia_ar5')
[7563]	88
	89	!Call to init moved to here so that we can call iom_use in the
	90	!initialisation
	91	IF( kt == nit000 ) CALL dia_ar5_init
[3294]	92
[11035]	93	CALL wrk_alloc( jpi , jpj , zarea_ssh , zbotpres, zpe)
[3294]	94	CALL wrk_alloc( jpi , jpj , jpk , zrhd , zrhop )
	95	CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn )
[1756]	96
	97	zarea_ssh(:,:) = area(:,:) * sshn(:,:)
	98
	99	! ! total volume of liquid seawater
	100	zvolssh = SUM( zarea_ssh(:,:) )
	101	IF( lk_mpp ) CALL mpp_sum( zvolssh )
	102	zvol = vol0 + zvolssh
	103
	104	CALL iom_put( 'voltot', zvol )
	105	CALL iom_put( 'sshtot', zvolssh / area_tot )
[7179]	106	CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
[1756]	107
[3294]	108	!
[7563]	109	IF( iom_use('sshthster')) THEN
	110	ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) ! thermosteric ssh
	111	ztsn(:,:,:,jp_sal) = sn0(:,:,:)
	112	CALL eos( ztsn, zrhd, fsdept_n(:,:,:) ) ! now in situ density using initial salinity
	113	!
	114	zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice
	115	DO jk = 1, jpkm1
	116	zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
	117	END DO
	118	IF( .NOT.lk_vvl ) THEN
	119	IF ( ln_isfcav ) THEN
	120	DO ji=1,jpi
	121	DO jj=1,jpj
	122	zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
	123	END DO
[5120]	124	END DO
[7563]	125	ELSE
	126	zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
	127	END IF
[5120]	128	END IF
[7563]	129	!
	130	zarho = SUM( area(:,:) * zbotpres(:,:) )
	131	IF( lk_mpp ) CALL mpp_sum( zarho )
	132	zssh_steric = - zarho / area_tot
	133	CALL iom_put( 'sshthster', zssh_steric )
[7179]	134	ENDIF
[1756]	135
	136	! ! steric sea surface height
[4313]	137	CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) ) ! now in situ and potential density
[2528]	138	zrhop(:,:,jpk) = 0._wp
[1756]	139	CALL iom_put( 'rhop', zrhop )
	140	!
[2528]	141	zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice
[1756]	142	DO jk = 1, jpkm1
	143	zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
	144	END DO
[4990]	145	IF( .NOT.lk_vvl ) THEN
[5120]	146	IF ( ln_isfcav ) THEN
	147	DO ji=1,jpi
	148	DO jj=1,jpj
	149	zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
	150	END DO
[4990]	151	END DO
[5120]	152	ELSE
	153	zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
	154	END IF
[4990]	155	END IF
[1756]	156	!
	157	zarho = SUM( area(:,:) * zbotpres(:,:) )
	158	IF( lk_mpp ) CALL mpp_sum( zarho )
	159	zssh_steric = - zarho / area_tot
	160	CALL iom_put( 'sshsteric', zssh_steric )
	161
	162	! ! ocean bottom pressure
[2528]	163	zztmp = rau0 * grav * 1.e-4_wp ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
[1756]	164	zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
	165	CALL iom_put( 'botpres', zbotpres )
	166
	167	! ! Mean density anomalie, temperature and salinity
[2528]	168	ztemp = 0._wp
	169	zsal = 0._wp
[1756]	170	DO jk = 1, jpkm1
	171	DO jj = 1, jpj
	172	DO ji = 1, jpi
	173	zztmp = area(ji,jj) * fse3t(ji,jj,jk)
[3294]	174	ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
	175	zsal = zsal + zztmp * tsn(ji,jj,jk,jp_sal)
[1756]	176	END DO
	177	END DO
	178	END DO
[2528]	179	IF( .NOT.lk_vvl ) THEN
[5120]	180	IF ( ln_isfcav ) THEN
	181	DO ji=1,jpi
	182	DO jj=1,jpj
	183	ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)
	184	zsal = zsal + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)
	185	END DO
[4990]	186	END DO
[5120]	187	ELSE
[5121]	188	ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
	189	zsal = zsal + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
[5120]	190	END IF
[1756]	191	ENDIF
	192	IF( lk_mpp ) THEN
	193	CALL mpp_sum( ztemp )
	194	CALL mpp_sum( zsal )
	195	END IF
	196	!
	197	zmass = rau0 * ( zarho + zvol ) ! total mass of liquid seawater
	198	ztemp = ztemp / zvol ! potential temperature in liquid seawater
	199	zsal = zsal / zvol ! Salinity of liquid seawater
	200	!
	201	CALL iom_put( 'masstot', zmass )
	202	CALL iom_put( 'temptot', ztemp )
	203	CALL iom_put( 'saltot' , zsal )
[7179]	204
	205	IF( iom_use( 'tnpeo' )) THEN
	206	! Work done against stratification by vertical mixing
	207	! Exclude points where rn2 is negative as convection kicks in here and
	208	! work is not being done against stratification
[11035]	209	zpe(:,:) = 0._wp
[7179]	210	IF( lk_zdfddm ) THEN
[11035]	211	DO jk = 2, jpk
	212	DO jj = 1, jpj
	213	DO ji = 1, jpi
	214	IF( rn2(ji,jj,jk) > 0._wp ) THEN
[7179]	215
[11035]	216	zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) &
	217	& / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )
	218	!
	219	zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
	220	zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
	221	!
	222
	223	zpe(ji, jj) = zpe(ji, jj) &
	224	- grav * ( avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) &
	225	- fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ))
	226
	227	ENDIF
[7179]	228	ENDDO
	229	ENDDO
	230	ENDDO
	231	ELSE
[11035]	232	DO jk=1,jpk
	233	DO ji=1,jpi
	234	DO jj=1,jpj
	235	zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
[7179]	236	ENDDO
	237	ENDDO
	238	ENDDO
	239	ENDIF
[11035]	240	CALL lbc_lnk(zpe, 'T', 1._wp)
	241	CALL iom_put( 'tnpeo', zpe )
[7179]	242	ENDIF
[1756]	243	!
[11035]	244	CALL wrk_dealloc( jpi , jpj , zarea_ssh , zbotpres, zpe )
[3294]	245	CALL wrk_dealloc( jpi , jpj , jpk , zrhd , zrhop )
	246	CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn )
[2715]	247	!
[3294]	248	IF( nn_timing == 1 ) CALL timing_stop('dia_ar5')
	249	!
[1756]	250	END SUBROUTINE dia_ar5
	251
	252
	253	SUBROUTINE dia_ar5_init
	254	!!----------------------------------------------------------------------
	255	!! * ROUTINE dia_ar5_init *
	256	!!
[2528]	257	!! ** Purpose : initialization for AR5 diagnostic computation
[1756]	258	!!----------------------------------------------------------------------
	259	INTEGER :: inum
	260	INTEGER :: ik
	261	INTEGER :: ji, jj, jk ! dummy loop indices
	262	REAL(wp) :: zztmp
[2715]	263	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zsaldta ! Jan/Dec levitus salinity
[5253]	264	!
[1756]	265	!!----------------------------------------------------------------------
	266	!
[3294]	267	IF( nn_timing == 1 ) CALL timing_start('dia_ar5_init')
	268	!
[6793]	269	CALL wrk_alloc( jpi, jpj, jpk, 2, zsaldta )
[2715]	270	! ! allocate dia_ar5 arrays
	271	IF( dia_ar5_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
	272
[1756]	273	area(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
	274
	275	area_tot = SUM( area(:,:) ) ; IF( lk_mpp ) CALL mpp_sum( area_tot )
	276
[2528]	277	vol0 = 0._wp
	278	thick0(:,:) = 0._wp
[1756]	279	DO jk = 1, jpkm1
[4292]	280	vol0 = vol0 + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) )
	281	thick0(:,:) = thick0(:,:) + tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk)
[1756]	282	END DO
[1948]	283	IF( lk_mpp ) CALL mpp_sum( vol0 )
[5253]	284
[7563]	285	IF( iom_use('sshthster')) THEN
	286	CALL iom_open ( 'sali_ref_clim_monthly', inum )
	287	CALL iom_get ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1 )
	288	CALL iom_get ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 )
	289	CALL iom_close( inum )
[6793]	290
[7563]	291	sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )
	292	sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
	293	IF( ln_zps ) THEN ! z-coord. partial steps
	294	DO jj = 1, jpj ! interpolation of salinity at the last ocean level (i.e. the partial step)
	295	DO ji = 1, jpi
	296	ik = mbkt(ji,jj)
	297	IF( ik > 1 ) THEN
	298	zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
	299	sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
	300	ENDIF
	301	END DO
[1756]	302	END DO
[7563]	303	ENDIF
[1756]	304	ENDIF
	305	!
[6793]	306	CALL wrk_dealloc( jpi, jpj, jpk, 2, zsaldta )
[2715]	307	!
[3294]	308	IF( nn_timing == 1 ) CALL timing_stop('dia_ar5_init')
	309	!
[1756]	310	END SUBROUTINE dia_ar5_init
	311
	312	#else
	313	!!----------------------------------------------------------------------
	314	!! Default option : NO diaar5
	315	!!----------------------------------------------------------------------
	316	LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .FALSE. ! coupled flag
	317	CONTAINS
[2528]	318	SUBROUTINE dia_ar5_init ! Dummy routine
	319	END SUBROUTINE dia_ar5_init
[1756]	320	SUBROUTINE dia_ar5( kt ) ! Empty routine
[2528]	321	INTEGER :: kt
[1756]	322	WRITE(,) 'dia_ar5: You should not have seen this print! error?', kt
	323	END SUBROUTINE dia_ar5
	324	#endif
	325
	326	!!======================================================================
	327	END MODULE diaar5

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/DIA/diaar5.F90

Download in other formats: