Context Navigation

iscplhsb.F90 @ 7753

Last change on this file since 7753 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

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

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

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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90 @ 7753

Download in other formats: