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

Last change on this file since 1321 was 1321, checked in by raasch, 11 years ago

last commit documented

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