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.

iscplhsb.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplhsb.F90 @ 10978

Last change on this file since 10978 was 10978, checked in by davestorkey, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : DOM and sshwzv.F90. (sshwzv.F90 will be rewritten somewhat when we rewrite the time-level swapping but I've done it anyway). Passes all non-AGRIF SETTE tests.
Property svn:keywords set to `Id`
File size: 16.1 KB

Rev	Line
[5790]	1	MODULE iscplhsb
	2	!!======================================================================
[7646]	3	!! * MODULE iscplhsb *
[5790]	4	!! Ocean forcing: ice sheet/ocean coupling (conservation)
	5	!!=====================================================================
	6	!! History : NEMO ! 2015-01 P. Mathiot: original
	7	!!----------------------------------------------------------------------
	8
	9	!!----------------------------------------------------------------------
	10	!! iscpl_alloc : variable allocation
	11	!! iscpl_hsb : compute and store the input of heat/salt/volume
	12	!! into the system due to the coupling process
	13	!! iscpl_div : correction of divergence to keep volume conservation
	14	!!----------------------------------------------------------------------
[9019]	15	USE oce ! global tra/dyn variable
[5790]	16	USE dom_oce ! ocean space and time domain
	17	USE domwri ! ocean space and time domain
[9019]	18	USE domngb !
[5790]	19	USE phycst ! physical constants
	20	USE sbc_oce ! surface boundary condition variables
[9019]	21	USE iscplini !
	22	!
[5790]	23	USE in_out_manager ! I/O manager
	24	USE lib_mpp ! MPP library
	25	USE lib_fortran ! MPP library
	26	USE lbclnk !
	27
	28	IMPLICIT NONE
	29	PRIVATE
	30
	31	PUBLIC iscpl_div
	32	PUBLIC iscpl_cons
	33	!! * Substitutions
[5920]	34	# include "vectopt_loop_substitute.h90"
[5790]	35	!!----------------------------------------------------------------------
[9598]	36	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10069]	37	!! $Id$
[10068]	38	!! Software governed by the CeCILL license (see ./LICENSE)
[5790]	39	!!----------------------------------------------------------------------
	40	CONTAINS
	41
[10978]	42	SUBROUTINE iscpl_cons( Kbb, Kmm, ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl )
[5790]	43	!!----------------------------------------------------------------------
	44	!! * ROUTINE iscpl_cons *
	45	!!
	46	!! ** Purpose : compute input into the system during the coupling step
	47	!! compute the correction term
	48	!! compute where the correction have to be applied
	49	!!
[10978]	50	!! ** Method : compute tsne3tn-tsbe3tb and e3tn-e3tb
[5790]	51	!!----------------------------------------------------------------------
[10978]	52	INTEGER , INTENT(in ) :: Kbb, Kmm !! time level indices
[5790]	53	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: ptmask_b !! mask before
	54	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pe3t_b !! scale factor before
	55	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: psmask_b !! mask before
	56	REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx !! corrective flux to have tracer conservation
	57	REAL(wp), DIMENSION(:,:,: ), INTENT(out) :: pvol_flx !! corrective flux to have volume conservation
	58	REAL(wp), INTENT(in ) :: prdt_iscpl !! coupling period
[9019]	59	!
	60	INTEGER :: ji , jj , jk ! loop index
	61	INTEGER :: jip1, jim1, jjp1, jjm1
	62	REAL(wp) :: summsk, zsum , zsumn, zjip1_ratio , zjim1_ratio, zdtem, zde3t, z1_rdtiscpl
	63	REAL(wp) :: zarea , zsum1, zsumb, zjjp1_ratio , zjjm1_ratio, zdsal
	64	REAL(wp), DIMENSION(jpi,jpj) :: zdssh ! workspace
	65	REAL(wp), DIMENSION(:), ALLOCATABLE :: zlon, zlat
	66	REAL(wp), DIMENSION(:), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
	67	INTEGER , DIMENSION(:), ALLOCATABLE :: ixpts, iypts, izpts, inpts
[5790]	68	INTEGER :: jpts, npts
[9019]	69	!!----------------------------------------------------------------------
[5790]	70
[5920]	71	! get imbalance (volume heat and salt)
	72	! initialisation difference
[9019]	73	zde3t = 0._wp ; zdsal = 0._wp ; zdtem = 0._wp
[5920]	74
	75	! initialisation correction term
[9019]	76	pvol_flx(:,:,: ) = 0._wp
	77	pts_flx (:,:,:,:) = 0._wp
[5920]	78
[9019]	79	z1_rdtiscpl = 1._wp / prdt_iscpl
[5790]	80
[10978]	81	! mask ts(:,:,:,:,Kmm) and ts(:,:,:,:,Kbb)
	82	ts(:,:,:,jp_tem,Kbb) = ts(:,:,:,jp_tem,Kbb) * ptmask_b(:,:,:)
	83	ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) * tmask (:,:,:)
	84	ts(:,:,:,jp_sal,Kbb) = ts(:,:,:,jp_sal,Kbb) * ptmask_b(:,:,:)
	85	ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) * tmask (:,:,:)
[5790]	86
[5835]	87	!==============================================================================
	88	! diagnose the heat, salt and volume input and compute the correction variable
	89	!==============================================================================
[5790]	90
[5835]	91	!
[10978]	92	zdssh(:,:) = ssh(:,:,Kmm) * ssmask(:,:) - ssh(:,:,Kbb) * psmask_b(:,:)
[6075]	93	IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
[5835]	94
[5802]	95	DO jk = 1,jpk-1
[5920]	96	DO jj = 2,jpj-1
	97	DO ji = fs_2,fs_jpim1
[5802]	98	IF (tmask_h(ji,jj) == 1._wp) THEN
[5790]	99
[5820]	100	! volume differences
[10978]	101	zde3t = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
[5790]	102
[5820]	103	! heat diff
[10978]	104	zdtem = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
	105	- ts(ji,jj,jk,jp_tem,Kbb) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
[5820]	106	! salt diff
[10978]	107	zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
	108	- ts(ji,jj,jk,jp_sal,Kbb) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
[5802]	109
[5820]	110	! shh changes
	111	IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN
[5920]	112	zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
	113	zdssh(ji,jj) = 0._wp
[5820]	114	END IF
[5790]	115
[5820]	116	! volume, heat and salt differences in each cell
[9019]	117	pvol_flx(ji,jj,jk) = pvol_flx(ji,jj,jk) + zde3t * z1_rdtiscpl
	118	pts_flx (ji,jj,jk,jp_sal)= pts_flx (ji,jj,jk,jp_sal) + zdsal * z1_rdtiscpl
	119	pts_flx (ji,jj,jk,jp_tem)= pts_flx (ji,jj,jk,jp_tem) + zdtem * z1_rdtiscpl
[5790]	120
[5835]	121	! case where we close a cell: check if the neighbour cells are wet
[5820]	122	IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
[5790]	123
[5820]	124	jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
[5790]	125
[5964]	126	zsum = e1e2t(ji ,jjp1) * tmask(ji ,jjp1,jk) + e1e2t(ji ,jjm1) * tmask(ji ,jjm1,jk) &
	127	& + e1e2t(jim1,jj ) * tmask(jim1,jj ,jk) + e1e2t(jip1,jj ) * tmask(jip1,jj ,jk)
[5790]	128
[5945]	129	IF ( zsum /= 0._wp ) THEN
[5964]	130	zjip1_ratio = e1e2t(jip1,jj ) * tmask(jip1,jj ,jk) / zsum
	131	zjim1_ratio = e1e2t(jim1,jj ) * tmask(jim1,jj ,jk) / zsum
	132	zjjp1_ratio = e1e2t(ji ,jjp1) * tmask(ji ,jjp1,jk) / zsum
	133	zjjm1_ratio = e1e2t(ji ,jjm1) * tmask(ji ,jjm1,jk) / zsum
[5790]	134
[5820]	135	pvol_flx(ji ,jjp1,jk ) = pvol_flx(ji ,jjp1,jk ) + pvol_flx(ji,jj,jk ) * zjjp1_ratio
	136	pvol_flx(ji ,jjm1,jk ) = pvol_flx(ji ,jjm1,jk ) + pvol_flx(ji,jj,jk ) * zjjm1_ratio
	137	pvol_flx(jip1,jj ,jk ) = pvol_flx(jip1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjip1_ratio
	138	pvol_flx(jim1,jj ,jk ) = pvol_flx(jim1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjim1_ratio
	139	pts_flx (ji ,jjp1,jk,jp_sal) = pts_flx (ji ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
	140	pts_flx (ji ,jjm1,jk,jp_sal) = pts_flx (ji ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
	141	pts_flx (jip1,jj ,jk,jp_sal) = pts_flx (jip1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
	142	pts_flx (jim1,jj ,jk,jp_sal) = pts_flx (jim1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
	143	pts_flx (ji ,jjp1,jk,jp_tem) = pts_flx (ji ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
	144	pts_flx (ji ,jjm1,jk,jp_tem) = pts_flx (ji ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
	145	pts_flx (jip1,jj ,jk,jp_tem) = pts_flx (jip1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
	146	pts_flx (jim1,jj ,jk,jp_tem) = pts_flx (jim1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
[5790]	147
[5802]	148	! set to 0 the cell we distributed over neigbourg cells
	149	pvol_flx(ji,jj,jk ) = 0._wp
	150	pts_flx (ji,jj,jk,jp_sal) = 0._wp
	151	pts_flx (ji,jj,jk,jp_tem) = 0._wp
[5820]	152
[5945]	153	ELSE IF (zsum == 0._wp ) THEN
[5820]	154	! case where we close a cell and no adjacent cell open
	155	! check if the cell beneath is wet
[5945]	156	IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
[5820]	157	pvol_flx(ji,jj,jk+1) = pvol_flx(ji,jj,jk+1) + pvol_flx(ji,jj,jk)
	158	pts_flx (ji,jj,jk+1,jp_sal)= pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
	159	pts_flx (ji,jj,jk+1,jp_tem)= pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
	160
	161	! set to 0 the cell we distributed over neigbourg cells
	162	pvol_flx(ji,jj,jk ) = 0._wp
	163	pts_flx (ji,jj,jk,jp_sal) = 0._wp
	164	pts_flx (ji,jj,jk,jp_tem) = 0._wp
	165	ELSE
	166	! case no adjacent cell on the horizontal and on the vertical
[5945]	167	IF ( lwp ) THEN ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
	168	! open for narea== 1 (lwp=T)
[5820]	169	WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
	170	WRITE(numout,*) ' ',mig(ji),' ',mjg(jj),' ',jk
	171	WRITE(numout,*) ' ',ji,' ',jj,' ',jk,' ',narea
	172	WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
[5945]	173	ENDIF
[5820]	174	! We deal with these points later.
	175	END IF
[5802]	176	END IF
	177	END IF
	178	END IF
	179	END DO
	180	END DO
	181	END DO
[5790]	182
[9019]	183	!!gm ERROR !!!!
	184	!! juste use tmask_i or in case of ISF smask_i (to be created to compute the sum without halos)
	185	!
	186	! CALL lbc_sum(pvol_flx(:,:,: ),'T',1.)
	187	! CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
	188	! CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
	189	STOP ' iscpl_cons: please modify this module !'
	190	!!gm end
[5802]	191	! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point
	192	! allocation and initialisation of the list of problematic point
[9019]	193	ALLOCATE( inpts(jpnij) )
	194	inpts(:) = 0
[5790]	195
[5802]	196	! fill narea location with the number of problematic point
	197	DO jk = 1,jpk-1
[5920]	198	DO jj = 2,jpj-1
	199	DO ji = fs_2,fs_jpim1
[5820]	200	IF ( ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1) == 0._wp .AND. tmask_h(ji,jj) == 1._wp &
	201	.AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
[5945]	202	inpts(narea) = inpts(narea) + 1
[5802]	203	END IF
	204	END DO
	205	END DO
	206	END DO
[5790]	207
[5802]	208	! build array of total problematic point on each cpu (share to each cpu)
[10425]	209	CALL mpp_max('iscplhsb', inpts,jpnij)
[5790]	210
[5802]	211	! size of the new variable
[5945]	212	npts = SUM(inpts)
[5802]	213
	214	! allocation of the coordinates, correction, index vector for the problematic points
	215	ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
[5945]	216	ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20_wp ; zlat(:) = -1.0e20_wp
	217	zcorr_vol(:) = -1.0e20_wp
	218	zcorr_sal(:) = -1.0e20_wp
	219	zcorr_tem(:) = -1.0e20_wp
[5790]	220
[5802]	221	! fill new variable
[5945]	222	jpts = SUM(inpts(1:narea-1))
[5802]	223	DO jk = 1,jpk-1
[5920]	224	DO jj = 2,jpj-1
	225	DO ji = fs_2,fs_jpim1
[5820]	226	IF ( ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1) == 0._wp .AND. tmask_h(ji,jj) == 1._wp &
	227	.AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
[5945]	228	jpts = jpts + 1 ! positioning in the inpts vector for the area narea
[5802]	229	ixpts(jpts) = ji ; iypts(jpts) = jj ; izpts(jpts) = jk
	230	zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
	231	zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
	232	zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
	233	zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
[5820]	234
[5802]	235	! set flx to 0 (safer)
	236	pvol_flx(ji,jj,jk ) = 0.0_wp
	237	pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
	238	pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
	239	END IF
	240	END DO
	241	END DO
	242	END DO
[5790]	243
[5802]	244	! build array of total problematic point on each cpu (share to each cpu)
[5820]	245	! point coordinates
[10425]	246	CALL mpp_max('iscplhsb', zlat ,npts)
	247	CALL mpp_max('iscplhsb', zlon ,npts)
	248	CALL mpp_max('iscplhsb', izpts,npts)
[5802]	249
[5820]	250	! correction values
[10425]	251	CALL mpp_max('iscplhsb', zcorr_vol,npts)
	252	CALL mpp_max('iscplhsb', zcorr_sal,npts)
	253	CALL mpp_max('iscplhsb', zcorr_tem,npts)
[5820]	254
[5802]	255	! put correction term in the closest cell
	256	DO jpts = 1,npts
	257	CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
	258	DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
	259	DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
	260	jk = izpts(jpts)
[5820]	261
	262	IF (tmask_h(ji,jj) == 1._wp) THEN
	263	! correct the vol_flx in the closest cell
	264	pvol_flx(ji,jj,jk) = pvol_flx(ji,jj,jk ) + zcorr_vol(jpts)
	265	pts_flx (ji,jj,jk,jp_sal) = pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
	266	pts_flx (ji,jj,jk,jp_tem) = pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
	267
	268	! set correction to 0
	269	zcorr_vol(jpts) = 0.0_wp
	270	zcorr_sal(jpts) = 0.0_wp
	271	zcorr_tem(jpts) = 0.0_wp
	272	END IF
[5802]	273	END DO
	274	END DO
	275	END DO
[5820]	276
[5802]	277	! deallocate variables
[5945]	278	DEALLOCATE(inpts)
[5802]	279	DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
	280
	281	! add contribution store on the hallo (lbclnk remove one of the contribution)
	282	pvol_flx(:,:,: ) = pvol_flx(:,:,: ) * tmask(:,:,:)
	283	pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
	284	pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
	285
[9019]	286	!!gm ERROR !!!!
	287	!! juste use tmask_i or in case of ISF smask_i (to be created to compute the sum without halos)
	288	!
	289	! ! compute sum over the halo and set it to 0.
	290	! CALL lbc_sum(pvol_flx(:,:,: ),'T',1._wp)
	291	! CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
	292	! CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
	293	!!gm end
[5802]	294
[9019]	295	!
[5790]	296	END SUBROUTINE iscpl_cons
	297
[9019]	298
[10978]	299	SUBROUTINE iscpl_div( Kmm, phdivn )
[5790]	300	!!----------------------------------------------------------------------
	301	!! * ROUTINE iscpl_div *
	302	!!
	303	!! ** Purpose : update the horizontal divergenc
	304	!!
	305	!! ** Method :
	306	!! CAUTION : iscpl is positive (inflow) and expressed in m/s
	307	!!
	308	!! ** Action : phdivn increase by the iscpl correction term
	309	!!----------------------------------------------------------------------
[10978]	310	INTEGER , INTENT(in ) :: Kmm ! time level index
[5790]	311	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence
	312	!!
	313	INTEGER :: ji, jj, jk ! dummy loop indices
	314	!!----------------------------------------------------------------------
	315	!
	316	DO jk = 1, jpk
	317	DO jj = 1, jpj
	318	DO ji = 1, jpi
[10978]	319	phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t(ji,jj,jk,Kmm)
[5790]	320	END DO
	321	END DO
	322	END DO
	323	!
	324	END SUBROUTINE iscpl_div
	325
	326	END MODULE iscplhsb

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

Download in other formats: