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

Last change on this file since 881 was 881, checked in by raasch, 9 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 raasch $
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(j*dny+jj,i*dnx+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(j*dny+jj,i*dnx+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.