Home

Context Navigation

source: palm/trunk/SOURCE/surface_coupler.f90 @ 1028

Last change on this file since 1028 was 881, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 18.7 KB

Line
1	SUBROUTINE surface_coupler
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! ------------------
10	! $Id: surface_coupler.f90 881 2012-04-13 06:36:58Z suehring $
11	!
12	! 880 2012-04-13 06:28:59Z raasch
13	! Bugfix: preprocessor statements for parallel execution added
14	!
15	! 709 2011-03-30 09:31:40Z raasch
16	! formatting adjustments
17	!
18	! 667 2010-12-23 12:06:00Z suehring/gryschka
19	! Additional case for nonequivalent processor and grid topopolgy in ocean and
20	! atmosphere added (coupling_topology = 1).
21	! Added exchange of u and v from Ocean to Atmosphere
22	!
23	! 291 2009-04-16 12:07:26Z raasch
24	! Coupling with independent precursor runs.
25	! Output of messages replaced by message handling routine.
26	!
27	! 206 2008-10-13 14:59:11Z raasch
28	! Implementation of a MPI-1 Coupling: replaced myid with target_id,
29	! deleted __mpi2 directives
30	!
31	! 109 2007-08-28 15:26:47Z letzel
32	! Initial revision
33	!
34	! Description:
35	! ------------
36	! Data exchange at the interface between coupled models
37	!------------------------------------------------------------------------------!
38
39	USE arrays_3d
40	USE control_parameters
41	USE cpulog
42	USE grid_variables
43	USE indices
44	USE interfaces
45	USE pegrid
46
47	IMPLICIT NONE
48
49	INTEGER :: i, j, k
50
51	REAL :: time_since_reference_point_rem
52	REAL :: total_2d(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
53
54	#if defined( __parallel )
55
56	CALL cpu_log( log_point(39), 'surface_coupler', 'start' )
57
58
59
60	!
61	!-- In case of model termination initiated by the remote model
62	!-- (terminate_coupled_remote > 0), initiate termination of the local model.
63	!-- The rest of the coupler must then be skipped because it would cause an MPI
64	!-- intercomminucation hang.
65	!-- If necessary, the coupler will be called at the beginning of the next
66	!-- restart run.
67
68	IF ( coupling_topology == 0 ) THEN
69	CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, target_id, &
70	0, &
71	terminate_coupled_remote, 1, MPI_INTEGER, target_id, &
72	0, comm_inter, status, ierr )
73	ELSE
74	IF ( myid == 0) THEN
75	CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, &
76	target_id, 0, &
77	terminate_coupled_remote, 1, MPI_INTEGER, &
78	target_id, 0, &
79	comm_inter, status, ierr )
80	ENDIF
81	CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
82	ierr )
83
84	ALLOCATE( total_2d_a(-nbgp:ny_a+nbgp,-nbgp:nx_a+nbgp), &
85	total_2d_o(-nbgp:ny_o+nbgp,-nbgp:nx_o+nbgp) )
86
87	ENDIF
88
89	IF ( terminate_coupled_remote > 0 ) THEN
90	WRITE( message_string, * ) 'remote model "', &
91	TRIM( coupling_mode_remote ), &
92	'" terminated', &
93	'&with terminate_coupled_remote = ', &
94	terminate_coupled_remote, &
95	'&local model "', TRIM( coupling_mode ), &
96	'" has', &
97	'&terminate_coupled = ', &
98	terminate_coupled
99	CALL message( 'surface_coupler', 'PA0310', 1, 2, 0, 6, 0 )
100	RETURN
101	ENDIF
102
103
104	!
105	!-- Exchange the current simulated time between the models,
106	!-- currently just for total_2ding
107	IF ( coupling_topology == 0 ) THEN
108
109	CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, 11, &
110	comm_inter, ierr )
111	CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, target_id, &
112	11, comm_inter, status, ierr )
113	ELSE
114
115	IF ( myid == 0 ) THEN
116
117	CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
118	11, comm_inter, ierr )
119	CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, &
120	target_id, 11, comm_inter, status, ierr )
121
122	ENDIF
123
124	CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
125	ierr )
126
127	ENDIF
128
129	!
130	!-- Exchange the interface data
131	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
132
133	!
134	!-- Horizontal grid size and number of processors is equal in ocean and
135	!-- atmosphere
136	IF ( coupling_topology == 0 ) THEN
137
138	!
139	!-- Send heat flux at bottom surface to the ocean
140	CALL MPI_SEND( shf(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
141	comm_inter, ierr )
142	!
143	!-- Send humidity flux at bottom surface to the ocean
144	IF ( humidity ) THEN
145	CALL MPI_SEND( qsws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 13, &
146	comm_inter, ierr )
147	ENDIF
148	!
149	!-- Receive temperature at the bottom surface from the ocean
150	CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, target_id, 14, &
151	comm_inter, status, ierr )
152	!
153	!-- Send the momentum flux (u) at bottom surface to the ocean
154	CALL MPI_SEND( usws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
155	comm_inter, ierr )
156	!
157	!-- Send the momentum flux (v) at bottom surface to the ocean
158	CALL MPI_SEND( vsws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
159	comm_inter, ierr )
160	!
161	!-- Receive u at the bottom surface from the ocean
162	CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, target_id, 17, &
163	comm_inter, status, ierr )
164	!
165	!-- Receive v at the bottom surface from the ocean
166	CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, target_id, 18, &
167	comm_inter, status, ierr )
168	!
169	!-- Horizontal grid size or number of processors differs between
170	!-- ocean and atmosphere
171	ELSE
172
173	!
174	!-- Send heat flux at bottom surface to the ocean
175	total_2d_a = 0.0
176	total_2d = 0.0
177	total_2d(nys:nyn,nxl:nxr) = shf(nys:nyn,nxl:nxr)
178
179	CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
180	comm2d, ierr )
181	CALL interpolate_to_ocean( 12 )
182	!
183	!-- Send humidity flux at bottom surface to the ocean
184	IF ( humidity ) THEN
185	total_2d_a = 0.0
186	total_2d = 0.0
187	total_2d(nys:nyn,nxl:nxr) = qsws(nys:nyn,nxl:nxr)
188
189	CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, &
190	0, comm2d, ierr )
191	CALL interpolate_to_ocean( 13 )
192	ENDIF
193	!
194	!-- Receive temperature at the bottom surface from the ocean
195	IF ( myid == 0 ) THEN
196	CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
197	target_id, 14, comm_inter, status, ierr )
198	ENDIF
199	CALL MPI_BARRIER( comm2d, ierr )
200	CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
201	ierr )
202	pt(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
203	!
204	!-- Send momentum flux (u) at bottom surface to the ocean
205	total_2d_a = 0.0
206	total_2d = 0.0
207	total_2d(nys:nyn,nxl:nxr) = usws(nys:nyn,nxl:nxr)
208	CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
209	comm2d, ierr )
210	CALL interpolate_to_ocean( 15 )
211	!
212	!-- Send momentum flux (v) at bottom surface to the ocean
213	total_2d_a = 0.0
214	total_2d = 0.0
215	total_2d(nys:nyn,nxl:nxr) = vsws(nys:nyn,nxl:nxr)
216	CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
217	comm2d, ierr )
218	CALL interpolate_to_ocean( 16 )
219	!
220	!-- Receive u at the bottom surface from the ocean
221	IF ( myid == 0 ) THEN
222	CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
223	target_id, 17, comm_inter, status, ierr )
224	ENDIF
225	CALL MPI_BARRIER( comm2d, ierr )
226	CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
227	ierr )
228	u(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
229	!
230	!-- Receive v at the bottom surface from the ocean
231	IF ( myid == 0 ) THEN
232	CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
233	target_id, 18, comm_inter, status, ierr )
234	ENDIF
235	CALL MPI_BARRIER( comm2d, ierr )
236	CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
237	ierr )
238	v(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
239
240	ENDIF
241
242	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
243
244	!
245	!-- Horizontal grid size and number of processors is equal
246	!-- in ocean and atmosphere
247	IF ( coupling_topology == 0 ) THEN
248	!
249	!-- Receive heat flux at the sea surface (top) from the atmosphere
250	CALL MPI_RECV( tswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
251	comm_inter, status, ierr )
252	!
253	!-- Receive humidity flux from the atmosphere (bottom)
254	!-- and add it to the heat flux at the sea surface (top)...
255	IF ( humidity_remote ) THEN
256	CALL MPI_RECV( qswst_remote(nysg,nxlg), ngp_xy, MPI_REAL, &
257	target_id, 13, comm_inter, status, ierr )
258	ENDIF
259	!
260	!-- Send sea surface temperature to the atmosphere model
261	CALL MPI_SEND( pt(nzt,nysg,nxlg), 1, type_xy, target_id, 14, &
262	comm_inter, ierr )
263	!
264	!-- Receive momentum flux (u) at the sea surface (top) from the atmosphere
265	CALL MPI_RECV( uswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
266	comm_inter, status, ierr )
267	!
268	!-- Receive momentum flux (v) at the sea surface (top) from the atmosphere
269	CALL MPI_RECV( vswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
270	comm_inter, status, ierr )
271	!
272	!-- Send u to the atmosphere
273	CALL MPI_SEND( u(nzt,nysg,nxlg), 1, type_xy, target_id, 17, &
274	comm_inter, ierr )
275	!
276	!-- Send v to the atmosphere
277	CALL MPI_SEND( v(nzt,nysg,nxlg), 1, type_xy, target_id, 18, &
278	comm_inter, ierr )
279	!
280	!-- Horizontal gridsize or number of processors differs between
281	!-- ocean and atmosphere
282	ELSE
283	!
284	!-- Receive heat flux at the sea surface (top) from the atmosphere
285	IF ( myid == 0 ) THEN
286	CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
287	target_id, 12, comm_inter, status, ierr )
288	ENDIF
289	CALL MPI_BARRIER( comm2d, ierr )
290	CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
291	ierr )
292	tswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
293	!
294	!-- Receive humidity flux at the sea surface (top) from the atmosphere
295	IF ( humidity_remote ) THEN
296	IF ( myid == 0 ) THEN
297	CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
298	target_id, 13, comm_inter, status, ierr )
299	ENDIF
300	CALL MPI_BARRIER( comm2d, ierr )
301	CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, &
302	comm2d, ierr)
303	qswst_remote(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
304	ENDIF
305	!
306	!-- Send surface temperature to atmosphere
307	total_2d_o = 0.0
308	total_2d = 0.0
309	total_2d(nys:nyn,nxl:nxr) = pt(nzt,nys:nyn,nxl:nxr)
310
311	CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
312	comm2d, ierr)
313	CALL interpolate_to_atmos( 14 )
314	!
315	!-- Receive momentum flux (u) at the sea surface (top) from the atmosphere
316	IF ( myid == 0 ) THEN
317	CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
318	target_id, 15, comm_inter, status, ierr )
319	ENDIF
320	CALL MPI_BARRIER( comm2d, ierr )
321	CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
322	0, comm2d, ierr )
323	uswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
324	!
325	!-- Receive momentum flux (v) at the sea surface (top) from the atmosphere
326	IF ( myid == 0 ) THEN
327	CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
328	target_id, 16, comm_inter, status, ierr )
329	ENDIF
330	CALL MPI_BARRIER( comm2d, ierr )
331	CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
332	ierr )
333	vswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg)
334	!
335	!-- Send u to atmosphere
336	total_2d_o = 0.0
337	total_2d = 0.0
338	total_2d(nys:nyn,nxl:nxr) = u(nzt,nys:nyn,nxl:nxr)
339	CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
340	comm2d, ierr )
341	CALL interpolate_to_atmos( 17 )
342	!
343	!-- Send v to atmosphere
344	total_2d_o = 0.0
345	total_2d = 0.0
346	total_2d(nys:nyn,nxl:nxr) = v(nzt,nys:nyn,nxl:nxr)
347	CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
348	comm2d, ierr )
349	CALL interpolate_to_atmos( 18 )
350
351	ENDIF
352
353	!
354	!-- Conversions of fluxes received from atmosphere
355	IF ( humidity_remote ) THEN
356	!
357	!-- Here tswst is still the sum of atmospheric bottom heat fluxes,
358	!-- * latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
359	!-- /(rho_atm(=1.0)*c_p)
360	tswst = tswst + qswst_remote * 2.2626108E6 / 1005.0
361	!
362	!-- ...and convert it to a salinity flux at the sea surface (top)
363	!-- following Steinhorn (1991), JPO 21, pp. 1681-1683:
364	!-- S'w' = -S * evaporation / ( rho_water * ( 1 - S ) )
365	saswst = -1.0 * sa(nzt,:,:) * qswst_remote / &
366	( rho(nzt,:,:) * ( 1.0 - sa(nzt,:,:) ) )
367	ENDIF
368
369	!
370	!-- Adjust the kinematic heat flux with respect to ocean density
371	!-- (constants are the specific heat capacities for air and water)
372	!-- now tswst is the ocean top heat flux
373	tswst = tswst / rho(nzt,:,:) * 1005.0 / 4218.0
374
375	!
376	!-- Adjust the momentum fluxes with respect to ocean density
377	uswst = uswst / rho(nzt,:,:)
378	vswst = vswst / rho(nzt,:,:)
379
380	ENDIF
381
382	IF ( coupling_topology == 1 ) THEN
383	DEALLOCATE( total_2d_o, total_2d_a )
384	ENDIF
385
386	CALL cpu_log( log_point(39), 'surface_coupler', 'stop' )
387
388	#endif
389
390	END SUBROUTINE surface_coupler
391
392
393
394	SUBROUTINE interpolate_to_atmos( tag )
395
396	#if defined( __parallel )
397
398	USE arrays_3d
399	USE control_parameters
400	USE grid_variables
401	USE indices
402	USE pegrid
403
404	IMPLICIT NONE
405
406	INTEGER :: dnx, dnx2, dny, dny2, i, ii, j, jj
407	INTEGER, intent(in) :: tag
408
409	CALL MPI_BARRIER( comm2d, ierr )
410
411	IF ( myid == 0 ) THEN
412	!
413	!-- Cyclic boundary conditions for the total 2D-grid
414	total_2d_o(-nbgp:-1,:) = total_2d_o(ny+1-nbgp:ny,:)
415	total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx+1-nbgp:nx)
416
417	total_2d_o(ny+1:ny+nbgp,:) = total_2d_o(0:nbgp-1,:)
418	total_2d_o(:,nx+1:nx+nbgp) = total_2d_o(:,0:nbgp-1)
419
420	!
421	!-- Number of gridpoints of the fine grid within one mesh of the coarse grid
422	dnx = (nx_o+1) / (nx_a+1)
423	dny = (ny_o+1) / (ny_a+1)
424
425	!
426	!-- Distance for interpolation around coarse grid points within the fine
427	!-- grid (note: 2*dnx2 must not be equal with dnx)
428	dnx2 = 2 * ( dnx / 2 )
429	dny2 = 2 * ( dny / 2 )
430
431	total_2d_a = 0.0
432	!
433	!-- Interpolation from ocean-grid-layer to atmosphere-grid-layer
434	DO j = 0, ny_a
435	DO i = 0, nx_a
436	DO jj = 0, dny2
437	DO ii = 0, dnx2
438	total_2d_a(j,i) = total_2d_a(j,i) &
439	+ total_2d_o(jdny+jj,idnx+ii)
440	ENDDO
441	ENDDO
442	total_2d_a(j,i) = total_2d_a(j,i) / ( ( dnx2 + 1 ) * ( dny2 + 1 ) )
443	ENDDO
444	ENDDO
445	!
446	!-- Cyclic boundary conditions for atmosphere grid
447	total_2d_a(-nbgp:-1,:) = total_2d_a(ny_a+1-nbgp:ny_a,:)
448	total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx_a+1-nbgp:nx_a)
449
450	total_2d_a(ny_a+1:ny_a+nbgp,:) = total_2d_a(0:nbgp-1,:)
451	total_2d_a(:,nx_a+1:nx_a+nbgp) = total_2d_a(:,0:nbgp-1)
452	!
453	!-- Transfer of the atmosphere-grid-layer to the atmosphere
454	CALL MPI_SEND( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, target_id, &
455	tag, comm_inter, ierr )
456
457	ENDIF
458
459	CALL MPI_BARRIER( comm2d, ierr )
460
461	#endif
462
463	END SUBROUTINE interpolate_to_atmos
464
465
466	SUBROUTINE interpolate_to_ocean( tag )
467
468	#if defined( __parallel )
469
470	USE arrays_3d
471	USE control_parameters
472	USE grid_variables
473	USE indices
474	USE pegrid
475
476	IMPLICIT NONE
477
478	INTEGER :: dnx, dny, i, ii, j, jj
479	INTEGER, intent(in) :: tag
480	REAL :: fl, fr, myl, myr
481
482
483	CALL MPI_BARRIER( comm2d, ierr )
484
485	IF ( myid == 0 ) THEN
486
487	!
488	!-- Number of gridpoints of the fine grid within one mesh of the coarse grid
489	dnx = ( nx_o + 1 ) / ( nx_a + 1 )
490	dny = ( ny_o + 1 ) / ( ny_a + 1 )
491
492	!
493	!-- Cyclic boundary conditions for atmosphere grid
494	total_2d_a(-nbgp:-1,:) = total_2d_a(ny+1-nbgp:ny,:)
495	total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx+1-nbgp:nx)
496
497	total_2d_a(ny+1:ny+nbgp,:) = total_2d_a(0:nbgp-1,:)
498	total_2d_a(:,nx+1:nx+nbgp) = total_2d_a(:,0:nbgp-1)
499	!
500	!-- Bilinear Interpolation from atmosphere grid-layer to ocean grid-layer
501	DO j = 0, ny
502	DO i = 0, nx
503	myl = ( total_2d_a(j+1,i) - total_2d_a(j,i) ) / dny
504	myr = ( total_2d_a(j+1,i+1) - total_2d_a(j,i+1) ) / dny
505	DO jj = 0, dny-1
506	fl = myl*jj + total_2d_a(j,i)
507	fr = myr*jj + total_2d_a(j,i+1)
508	DO ii = 0, dnx-1
509	total_2d_o(jdny+jj,idnx+ii) = ( fr - fl ) / dnx * ii + fl
510	ENDDO
511	ENDDO
512	ENDDO
513	ENDDO
514	!
515	!-- Cyclic boundary conditions for ocean grid
516	total_2d_o(-nbgp:-1,:) = total_2d_o(ny_o+1-nbgp:ny_o,:)
517	total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx_o+1-nbgp:nx_o)
518
519	total_2d_o(ny_o+1:ny_o+nbgp,:) = total_2d_o(0:nbgp-1,:)
520	total_2d_o(:,nx_o+1:nx_o+nbgp) = total_2d_o(:,0:nbgp-1)
521
522	CALL MPI_SEND( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
523	target_id, tag, comm_inter, ierr )
524
525	ENDIF
526
527	CALL MPI_BARRIER( comm2d, ierr )
528
529	#endif
530
531	END SUBROUTINE interpolate_to_ocean

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |