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

Last change on this file since 2550 was 2233, checked in by suehring, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 31.6 KB
Line 
1!> @file surface_coupler.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: surface_coupler.f90 2233 2017-05-30 18:08:54Z boeske $
27!
28! 2232 2017-05-30 17:47:52Z suehring
29! Adjust to new surface structure. Transfer 1D surface fluxes onto 2D grid
30! (and back).
31!
32! 2031 2016-10-21 15:11:58Z knoop
33! renamed variable rho to rho_ocean
34!
35! 2000 2016-08-20 18:09:15Z knoop
36! Forced header and separation lines into 80 columns
37!
38! 1682 2015-10-07 23:56:08Z knoop
39! Code annotations made doxygen readable
40!
41! 1427 2014-07-07 14:04:59Z maronga
42! Bugfix: value of l_v corrected.
43!
44! 1418 2014-06-06 13:05:08Z fricke
45! Bugfix: For caluclation of the salinity flux at the sea surface,
46!          the given value for salinity must be in percent and not in psu
47!
48! 1353 2014-04-08 15:21:23Z heinze
49! REAL constants provided with KIND-attribute
50!
51! 1324 2014-03-21 09:13:16Z suehring
52! Bugfix: ONLY statement for module pegrid removed
53!
54! 1322 2014-03-20 16:38:49Z raasch
55! REAL constants defined as wp-kind
56!
57! 1320 2014-03-20 08:40:49Z raasch
58! ONLY-attribute added to USE-statements,
59! kind-parameters added to all INTEGER and REAL declaration statements,
60! kinds are defined in new module kinds,
61! old module precision_kind is removed,
62! revision history before 2012 removed,
63! comment fields (!:) to be used for variable explanations added to
64! all variable declaration statements
65!
66! 1318 2014-03-17 13:35:16Z raasch
67! module interfaces removed
68!
69! 1092 2013-02-02 11:24:22Z raasch
70! unused variables removed
71!
72! 1036 2012-10-22 13:43:42Z raasch
73! code put under GPL (PALM 3.9)
74!
75! 880 2012-04-13 06:28:59Z raasch
76! Bugfix: preprocessor statements for parallel execution added
77!
78! 109 2007-08-28 15:26:47Z letzel
79! Initial revision
80!
81! Description:
82! ------------
83!> Data exchange at the interface between coupled models
84!------------------------------------------------------------------------------!
85 SUBROUTINE surface_coupler
86 
87
88    USE arrays_3d,                                                             &
89        ONLY:  pt, rho_ocean, sa, total_2d_a, total_2d_o, u, v
90
91    USE cloud_parameters,                                                      &
92        ONLY:  cp, l_v
93
94    USE control_parameters,                                                    &
95        ONLY:  coupling_mode, coupling_mode_remote, coupling_topology,         &
96               humidity, humidity_remote, land_surface, message_string,        &
97               terminate_coupled, terminate_coupled_remote,                    &
98               time_since_reference_point, urban_surface
99
100    USE cpulog,                                                                &
101        ONLY:  cpu_log, log_point
102
103    USE indices,                                                               &
104        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nx_a, nx_o, ny, nyn, nyng, nys, &
105               nysg, ny_a, ny_o, nzt
106
107    USE kinds
108
109    USE pegrid
110
111    USE surface_mod,                                                           &
112        ONLY :  surf_def_h, surf_lsm_h, surf_type, surf_usm_h
113
114    IMPLICIT NONE
115
116    INTEGER(iwp) ::  i                                    !< index variable x-direction
117    INTEGER(iwp) ::  j                                    !< index variable y-direction
118    INTEGER(iwp) ::  m                                    !< running index for surface elements
119
120    REAL(wp)    ::  cpw = 4218.0_wp                       !< heat capacity of water at constant pressure
121    REAL(wp)    ::  time_since_reference_point_rem        !<
122    REAL(wp)    ::  total_2d(-nbgp:ny+nbgp,-nbgp:nx+nbgp) !<
123
124    REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  surface_flux !< dummy array for surface fluxes on 2D grid
125
126
127#if defined( __parallel )
128
129    CALL cpu_log( log_point(39), 'surface_coupler', 'start' )
130
131
132
133!
134!-- In case of model termination initiated by the remote model
135!-- (terminate_coupled_remote > 0), initiate termination of the local model.
136!-- The rest of the coupler must then be skipped because it would cause an MPI
137!-- intercomminucation hang.
138!-- If necessary, the coupler will be called at the beginning of the next
139!-- restart run.
140
141    IF ( coupling_topology == 0 ) THEN
142       CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, target_id, &
143                          0,                                                   &
144                          terminate_coupled_remote, 1, MPI_INTEGER, target_id, &
145                          0, comm_inter, status, ierr )
146    ELSE
147       IF ( myid == 0) THEN
148          CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, &
149                             target_id, 0,                             &
150                             terminate_coupled_remote, 1, MPI_INTEGER, & 
151                             target_id, 0,                             &
152                             comm_inter, status, ierr )
153       ENDIF
154       CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
155                       ierr )
156
157       ALLOCATE( total_2d_a(-nbgp:ny_a+nbgp,-nbgp:nx_a+nbgp),       &
158                 total_2d_o(-nbgp:ny_o+nbgp,-nbgp:nx_o+nbgp) )
159
160    ENDIF
161
162    IF ( terminate_coupled_remote > 0 )  THEN
163       WRITE( message_string, * ) 'remote model "',                         &
164                                  TRIM( coupling_mode_remote ),             &
165                                  '" terminated',                           &
166                                  '&with terminate_coupled_remote = ',      &
167                                  terminate_coupled_remote,                 &
168                                  '&local model  "', TRIM( coupling_mode ), &
169                                  '" has',                                  &
170                                  '&terminate_coupled = ',                  &
171                                   terminate_coupled
172       CALL message( 'surface_coupler', 'PA0310', 1, 2, 0, 6, 0 )
173       RETURN
174    ENDIF
175 
176
177!
178!-- Exchange the current simulated time between the models,
179!-- currently just for total_2d
180    IF ( coupling_topology == 0 ) THEN
181   
182       CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, 11, &
183                      comm_inter, ierr )
184       CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, target_id, &
185                      11, comm_inter, status, ierr )
186    ELSE
187
188       IF ( myid == 0 ) THEN
189
190          CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
191                         11, comm_inter, ierr )
192          CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL,        &
193                         target_id, 11, comm_inter, status, ierr )
194
195       ENDIF
196
197       CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
198                       ierr )
199
200    ENDIF
201
202!
203!-- Exchange the interface data
204    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
205   
206!
207!--    Horizontal grid size and number of processors is equal in ocean and
208!--    atmosphere
209       IF ( coupling_topology == 0 )  THEN
210
211!
212!--       Send heat flux at bottom surface to the ocean. First, transfer from
213!--       1D surface type to 2D grid.
214          CALL transfer_1D_to_2D_equal( surf_def_h(0)%shf, surf_lsm_h%shf,     &
215                                        surf_usm_h%shf )
216          CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, &
217                         12, comm_inter, ierr )
218!
219!--       Send humidity flux at bottom surface to the ocean. First, transfer
220!--       from 1D surface type to 2D grid.
221          CALL transfer_1D_to_2D_equal( surf_def_h(0)%qsws, surf_lsm_h%qsws,   &
222                                        surf_usm_h%qsws )
223          IF ( humidity )  THEN
224             CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL,         &
225                            target_id, 13, comm_inter, ierr )
226          ENDIF
227!
228!--       Receive temperature at the bottom surface from the ocean
229          CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, target_id, 14,           &
230                         comm_inter, status, ierr )
231!
232!--       Send the momentum flux (u) at bottom surface to the ocean. First,
233!--       transfer from 1D surface type to 2D grid.
234          CALL transfer_1D_to_2D_equal( surf_def_h(0)%usws, surf_lsm_h%usws,   &
235                                        surf_usm_h%usws )
236          CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, &
237                         15, comm_inter, ierr )
238!
239!--       Send the momentum flux (v) at bottom surface to the ocean. First,
240!--       transfer from 1D surface type to 2D grid.
241          CALL transfer_1D_to_2D_equal( surf_def_h(0)%vsws, surf_lsm_h%vsws,   &
242                                        surf_usm_h%vsws )
243          CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, &
244                         16, comm_inter, ierr )
245!
246!--       Receive u at the bottom surface from the ocean
247          CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, target_id, 17,            &
248                         comm_inter, status, ierr )
249!
250!--       Receive v at the bottom surface from the ocean
251          CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, target_id, 18,            &
252                         comm_inter, status, ierr )
253!
254!--    Horizontal grid size or number of processors differs between
255!--    ocean and atmosphere
256       ELSE
257     
258!
259!--       Send heat flux at bottom surface to the ocean
260          total_2d_a = 0.0_wp
261          total_2d   = 0.0_wp
262!
263!--       Transfer from 1D surface type to 2D grid.
264          CALL transfer_1D_to_2D_unequal( surf_def_h(0)%shf, surf_lsm_h%shf,   &
265                                          surf_usm_h%shf )
266
267          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0,  &
268                           comm2d, ierr )
269          CALL interpolate_to_ocean( 12 )   
270!
271!--       Send humidity flux at bottom surface to the ocean
272          IF ( humidity )  THEN
273             total_2d_a = 0.0_wp
274             total_2d   = 0.0_wp
275!
276!--          Transfer from 1D surface type to 2D grid.
277             CALL transfer_1D_to_2D_unequal( surf_def_h(0)%qsws,              &
278                                             surf_lsm_h%qsws,                 &
279                                             surf_usm_h%qsws )
280
281             CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, &
282                              0, comm2d, ierr )
283             CALL interpolate_to_ocean( 13 )
284          ENDIF
285!
286!--       Receive temperature at the bottom surface from the ocean
287          IF ( myid == 0 )  THEN
288             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL,          &
289                            target_id, 14, comm_inter, status, ierr )   
290          ENDIF
291          CALL MPI_BARRIER( comm2d, ierr )
292          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
293                          ierr )
294          pt(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
295!
296!--       Send momentum flux (u) at bottom surface to the ocean
297          total_2d_a = 0.0_wp 
298          total_2d   = 0.0_wp
299!
300!--       Transfer from 1D surface type to 2D grid.
301          CALL transfer_1D_to_2D_unequal( surf_def_h(0)%usws, surf_lsm_h%usws, &
302                                          surf_usm_h%usws )
303          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
304                           comm2d, ierr )
305          CALL interpolate_to_ocean( 15 )
306!
307!--       Send momentum flux (v) at bottom surface to the ocean
308          total_2d_a = 0.0_wp
309          total_2d   = 0.0_wp
310!
311!--       Transfer from 1D surface type to 2D grid.
312          CALL transfer_1D_to_2D_unequal( surf_def_h(0)%usws, surf_lsm_h%usws, &
313                                          surf_usm_h%usws )
314          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
315                           comm2d, ierr )
316          CALL interpolate_to_ocean( 16 )
317!
318!--       Receive u at the bottom surface from the ocean
319          IF ( myid == 0 )  THEN
320             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
321                            target_id, 17, comm_inter, status, ierr )
322          ENDIF
323          CALL MPI_BARRIER( comm2d, ierr )
324          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
325                          ierr )
326          u(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
327!
328!--       Receive v at the bottom surface from the ocean
329          IF ( myid == 0 )  THEN
330             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
331                            target_id, 18, comm_inter, status, ierr )
332          ENDIF
333          CALL MPI_BARRIER( comm2d, ierr )
334          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
335                          ierr )
336          v(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
337
338       ENDIF
339
340    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
341
342!
343!--    Horizontal grid size and number of processors is equal
344!--    in ocean and atmosphere
345       IF ( coupling_topology == 0 ) THEN
346!
347!--       Receive heat flux at the sea surface (top) from the atmosphere
348          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
349                         comm_inter, status, ierr )
350          CALL transfer_2D_to_1D_equal( surf_def_h(2)%shf )
351!
352!--       Receive humidity flux from the atmosphere (bottom)
353!--       and add it to the heat flux at the sea surface (top)...
354          IF ( humidity_remote )  THEN
355             CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, &
356                            target_id, 13, comm_inter, status, ierr )
357             CALL transfer_2D_to_1D_equal( surf_def_h(2)%qsws )
358          ENDIF
359!
360!--       Send sea surface temperature to the atmosphere model
361          CALL MPI_SEND( pt(nzt,nysg,nxlg), 1, type_xy, target_id, 14, &
362                         comm_inter, ierr )
363!
364!--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
365          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
366                         comm_inter, status, ierr )
367          CALL transfer_2D_to_1D_equal( surf_def_h(2)%usws )
368!
369!--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
370          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
371                         comm_inter, status, ierr )
372          CALL transfer_2D_to_1D_equal( surf_def_h(2)%vsws )
373!
374!--       Send u to the atmosphere
375          CALL MPI_SEND( u(nzt,nysg,nxlg), 1, type_xy, target_id, 17, &
376                         comm_inter, ierr )
377!
378!--       Send v to the atmosphere
379          CALL MPI_SEND( v(nzt,nysg,nxlg), 1, type_xy, target_id, 18, &
380                         comm_inter, ierr )
381!
382!--    Horizontal gridsize or number of processors differs between
383!--    ocean and atmosphere
384       ELSE
385!
386!--       Receive heat flux at the sea surface (top) from the atmosphere
387          IF ( myid == 0 )  THEN
388             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
389                            target_id, 12, comm_inter, status, ierr )
390          ENDIF
391          CALL MPI_BARRIER( comm2d, ierr )
392          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
393                          ierr )
394          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%shf )
395!
396!--       Receive humidity flux at the sea surface (top) from the atmosphere
397          IF ( humidity_remote )  THEN
398             IF ( myid == 0 )  THEN
399                CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
400                               target_id, 13, comm_inter, status, ierr )
401             ENDIF
402             CALL MPI_BARRIER( comm2d, ierr )
403             CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, &
404                             comm2d, ierr)
405             CALL transfer_2D_to_1D_unequal( surf_def_h(2)%qsws )
406          ENDIF
407!
408!--       Send surface temperature to atmosphere
409          total_2d_o = 0.0_wp
410          total_2d   = 0.0_wp
411          total_2d(nys:nyn,nxl:nxr) = pt(nzt,nys:nyn,nxl:nxr)
412
413          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
414                           comm2d, ierr) 
415          CALL interpolate_to_atmos( 14 )
416!
417!--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
418          IF ( myid == 0 )  THEN
419             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
420                            target_id, 15, comm_inter, status, ierr )
421          ENDIF
422          CALL MPI_BARRIER( comm2d, ierr )
423          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
424                          0, comm2d, ierr )
425          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%usws )
426!
427!--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
428          IF ( myid == 0 )  THEN
429             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
430                            target_id, 16, comm_inter, status, ierr )
431          ENDIF
432          CALL MPI_BARRIER( comm2d, ierr )
433          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
434                          ierr )
435          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%vsws )
436!
437!--       Send u to atmosphere
438          total_2d_o = 0.0_wp 
439          total_2d   = 0.0_wp
440          total_2d(nys:nyn,nxl:nxr) = u(nzt,nys:nyn,nxl:nxr)
441          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
442                           comm2d, ierr )
443          CALL interpolate_to_atmos( 17 )
444!
445!--       Send v to atmosphere
446          total_2d_o = 0.0_wp
447          total_2d   = 0.0_wp
448          total_2d(nys:nyn,nxl:nxr) = v(nzt,nys:nyn,nxl:nxr)
449          CALL MPI_REDUCE( total_2d, total_2d_o, ngp_o, MPI_REAL, MPI_SUM, 0, &
450                           comm2d, ierr )
451          CALL interpolate_to_atmos( 18 )
452       
453       ENDIF
454
455!
456!--    Conversions of fluxes received from atmosphere
457       IF ( humidity_remote )  THEN
458!
459!--       Here top heat flux is still the sum of atmospheric bottom heat fluxes,
460!--       * latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
461!--       /(rho_atm(=1.0)*c_p)
462          DO  m = 1, surf_def_h(2)%ns
463             i = surf_def_h(2)%i(m)
464             j = surf_def_h(2)%j(m)
465             
466             surf_def_h(2)%shf(m) = surf_def_h(2)%shf(m) +                     &
467                                    surf_def_h(2)%qsws(m) * l_v / cp
468!
469!--          ...and convert it to a salinity flux at the sea surface (top)
470!--          following Steinhorn (1991), JPO 21, pp. 1681-1683:
471!--          S'w' = -S * evaporation / ( rho_water * ( 1 - S ) )
472             surf_def_h(2)%sasws(m) = -1.0_wp * sa(nzt,j,i) * 0.001_wp *       &
473                                      surf_def_h(2)%qsws(m) /                  &
474                                    ( rho_ocean(nzt,j,i) *                     &
475                                      ( 1.0_wp - sa(nzt,j,i) * 0.001_wp )      &
476                                    )
477          ENDDO
478       ENDIF
479
480!
481!--    Adjust the kinematic heat flux with respect to ocean density
482!--    (constants are the specific heat capacities for air and water), as well
483!--    as momentum fluxes
484       DO  m = 1, surf_def_h(2)%ns
485          i = surf_def_h(2)%i(m)
486          j = surf_def_h(2)%j(m)
487          surf_def_h(2)%shf(m) = surf_def_h(2)%shf(m) / rho_ocean(nzt,j,i) *   &
488                                 cp / cpw
489
490          surf_def_h(2)%usws(m) = surf_def_h(2)%usws(m) / rho_ocean(nzt,j,i)
491          surf_def_h(2)%vsws(m) = surf_def_h(2)%vsws(m) / rho_ocean(nzt,j,i)
492       ENDDO
493
494    ENDIF
495
496    IF ( coupling_topology == 1 )  THEN
497       DEALLOCATE( total_2d_o, total_2d_a )
498    ENDIF
499
500    CALL cpu_log( log_point(39), 'surface_coupler', 'stop' )
501
502#endif
503
504     CONTAINS 
505
506!       Description:
507!------------------------------------------------------------------------------!
508!>      Data transfer from 1D surface-data type to 2D dummy array for equal
509!>      grids in atmosphere and ocean.
510!------------------------------------------------------------------------------!
511        SUBROUTINE transfer_1D_to_2D_equal( def_1d, lsm_1d, usm_1d )
512
513           IMPLICIT NONE
514
515            INTEGER(iwp) ::  i   !< running index x
516            INTEGER(iwp) ::  j   !< running index y
517            INTEGER(iwp) ::  m   !< running index surface type
518
519            REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  def_1d !< 1D surface flux, default surfaces
520            REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  lsm_1d !< 1D surface flux, natural surfaces
521            REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  usm_1d !< 1D surface flux, urban surfaces
522!
523!--         Transfer surface flux at default surfaces to 2D grid
524            DO  m = 1, surf_def_h(0)%ns
525               i = surf_def_h(0)%i(m)
526               j = surf_def_h(0)%j(m)
527               surface_flux(j,i) = def_1d(m)
528            ENDDO
529!
530!--         Transfer surface flux at natural surfaces to 2D grid
531            IF ( land_surface )  THEN
532               DO  m = 1, SIZE(lsm_1d)
533                  i = surf_lsm_h%i(m)
534                  j = surf_lsm_h%j(m)
535                  surface_flux(j,i) = lsm_1d(m)
536               ENDDO
537            ENDIF
538!
539!--         Transfer surface flux at natural surfaces to 2D grid
540            IF ( urban_surface )  THEN
541               DO  m = 1, SIZE(usm_1d)
542                  i = surf_usm_h%i(m)
543                  j = surf_usm_h%j(m)
544                  surface_flux(j,i) = usm_1d(m)
545               ENDDO
546            ENDIF
547
548        END SUBROUTINE transfer_1D_to_2D_equal
549
550!       Description:
551!------------------------------------------------------------------------------!
552!>      Data transfer from 2D array for equal grids onto 1D surface-data type
553!>      array.
554!------------------------------------------------------------------------------!
555        SUBROUTINE transfer_2D_to_1D_equal( def_1d )
556
557           IMPLICIT NONE
558
559            INTEGER(iwp) ::  i   !< running index x
560            INTEGER(iwp) ::  j   !< running index y
561            INTEGER(iwp) ::  m   !< running index surface type
562
563            REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  def_1d !< 1D surface flux, default surfaces
564!
565!--         Transfer surface flux to 1D surface type, only for default surfaces
566            DO  m = 1, surf_def_h(2)%ns
567               i = surf_def_h(2)%i(m)
568               j = surf_def_h(2)%j(m)
569               def_1d(m) = surface_flux(j,i)
570            ENDDO
571
572        END SUBROUTINE transfer_2D_to_1D_equal
573
574!       Description:
575!------------------------------------------------------------------------------!
576!>      Data transfer from 1D surface-data type to 2D dummy array from unequal
577!>      grids in atmosphere and ocean.
578!------------------------------------------------------------------------------!
579        SUBROUTINE transfer_1D_to_2D_unequal( def_1d, lsm_1d, usm_1d )
580
581           IMPLICIT NONE
582
583            INTEGER(iwp) ::  i   !< running index x
584            INTEGER(iwp) ::  j   !< running index y
585            INTEGER(iwp) ::  m   !< running index surface type
586
587            REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  def_1d !< 1D surface flux, default surfaces
588            REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  lsm_1d !< 1D surface flux, natural surfaces
589            REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  usm_1d !< 1D surface flux, urban surfaces
590!
591!--         Transfer surface flux at default surfaces to 2D grid. Transfer no
592!--         ghost-grid points since total_2d is a global array.
593            DO  m = 1, SIZE(def_1d)
594               i = surf_def_h(0)%i(m)
595               j = surf_def_h(0)%j(m)
596
597               IF ( i >= nxl  .AND.  i <= nxr  .AND.                           &
598                    j >= nys  .AND.  j <= nyn )  THEN
599                  total_2d(j,i) = def_1d(m)
600               ENDIF
601            ENDDO
602!
603!--         Transfer surface flux at natural surfaces to 2D grid
604            IF ( land_surface )  THEN
605               DO  m = 1, SIZE(lsm_1d)
606                  i = surf_lsm_h%i(m)
607                  j = surf_lsm_h%j(m)
608
609                  IF ( i >= nxl  .AND.  i <= nxr  .AND.                        &
610                       j >= nys  .AND.  j <= nyn )  THEN
611                     total_2d(j,i) = lsm_1d(m)
612                  ENDIF
613               ENDDO
614            ENDIF
615!
616!--         Transfer surface flux at natural surfaces to 2D grid
617            IF ( urban_surface )  THEN
618               DO  m = 1, SIZE(usm_1d)
619                  i = surf_usm_h%i(m)
620                  j = surf_usm_h%j(m)
621
622                  IF ( i >= nxl  .AND.  i <= nxr  .AND.                        &
623                       j >= nys  .AND.  j <= nyn )  THEN
624                     total_2d(j,i) = usm_1d(m)
625                  ENDIF
626               ENDDO
627            ENDIF
628
629        END SUBROUTINE transfer_1D_to_2D_unequal
630
631!       Description:
632!------------------------------------------------------------------------------!
633!>      Data transfer from 2D dummy array from unequal grids to 1D surface-data
634!>      type.
635!------------------------------------------------------------------------------!
636        SUBROUTINE transfer_2D_to_1D_unequal( def_1d )
637
638           IMPLICIT NONE
639
640            INTEGER(iwp) ::  i   !< running index x
641            INTEGER(iwp) ::  j   !< running index y
642            INTEGER(iwp) ::  m   !< running index surface type
643
644            REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  def_1d !< 1D surface flux, default surfaces
645!
646!--         Transfer 2D surface flux to default surfaces data type. Transfer no
647!--         ghost-grid points since total_2d is a global array.
648            DO  m = 1, SIZE(def_1d)
649               i = surf_def_h(2)%i(m)
650               j = surf_def_h(2)%j(m)
651
652               IF ( i >= nxl  .AND.  i <= nxr  .AND.                           &
653                    j >= nys  .AND.  j <= nyn )  THEN
654                  def_1d(m) = total_2d_o(j,i)
655               ENDIF
656            ENDDO
657
658
659        END SUBROUTINE transfer_2D_to_1D_unequal
660
661  END SUBROUTINE surface_coupler
662
663
664
665!------------------------------------------------------------------------------!
666! Description:
667! ------------
668!> @todo Missing subroutine description.
669!------------------------------------------------------------------------------!
670  SUBROUTINE interpolate_to_atmos( tag )
671
672#if defined( __parallel )
673
674    USE arrays_3d,                                                             &
675        ONLY:  total_2d_a, total_2d_o
676
677    USE indices,                                                               &
678        ONLY:  nbgp, nx, nx_a, nx_o, ny, ny_a, ny_o
679
680    USE kinds
681
682    USE pegrid
683
684    IMPLICIT NONE
685
686    INTEGER(iwp) ::  dnx  !<
687    INTEGER(iwp) ::  dnx2 !<
688    INTEGER(iwp) ::  dny  !<
689    INTEGER(iwp) ::  dny2 !<
690    INTEGER(iwp) ::  i    !<
691    INTEGER(iwp) ::  ii   !<
692    INTEGER(iwp) ::  j    !<
693    INTEGER(iwp) ::  jj   !<
694
695    INTEGER(iwp), intent(in) ::  tag !<
696
697    CALL MPI_BARRIER( comm2d, ierr )
698
699    IF ( myid == 0 )  THEN
700!
701!--    Cyclic boundary conditions for the total 2D-grid
702       total_2d_o(-nbgp:-1,:) = total_2d_o(ny+1-nbgp:ny,:)
703       total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx+1-nbgp:nx)
704
705       total_2d_o(ny+1:ny+nbgp,:) = total_2d_o(0:nbgp-1,:)
706       total_2d_o(:,nx+1:nx+nbgp) = total_2d_o(:,0:nbgp-1)
707
708!
709!--    Number of gridpoints of the fine grid within one mesh of the coarse grid
710       dnx = (nx_o+1) / (nx_a+1) 
711       dny = (ny_o+1) / (ny_a+1) 
712
713!
714!--    Distance for interpolation around coarse grid points within the fine
715!--    grid (note: 2*dnx2 must not be equal with dnx)
716       dnx2 = 2 * ( dnx / 2 )
717       dny2 = 2 * ( dny / 2 )
718
719       total_2d_a = 0.0_wp
720!
721!--    Interpolation from ocean-grid-layer to atmosphere-grid-layer
722       DO  j = 0, ny_a
723          DO  i = 0, nx_a 
724             DO  jj = 0, dny2
725                DO  ii = 0, dnx2
726                   total_2d_a(j,i) = total_2d_a(j,i) &
727                                     + total_2d_o(j*dny+jj,i*dnx+ii)
728                ENDDO
729             ENDDO
730             total_2d_a(j,i) = total_2d_a(j,i) / ( ( dnx2 + 1 ) * ( dny2 + 1 ) )
731          ENDDO
732       ENDDO
733!
734!--    Cyclic boundary conditions for atmosphere grid
735       total_2d_a(-nbgp:-1,:) = total_2d_a(ny_a+1-nbgp:ny_a,:)
736       total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx_a+1-nbgp:nx_a)
737       
738       total_2d_a(ny_a+1:ny_a+nbgp,:) = total_2d_a(0:nbgp-1,:)
739       total_2d_a(:,nx_a+1:nx_a+nbgp) = total_2d_a(:,0:nbgp-1)
740!
741!--    Transfer of the atmosphere-grid-layer to the atmosphere
742       CALL MPI_SEND( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, target_id, &
743                      tag, comm_inter, ierr )
744
745    ENDIF
746
747    CALL MPI_BARRIER( comm2d, ierr )
748
749#endif
750
751  END SUBROUTINE interpolate_to_atmos
752
753
754!------------------------------------------------------------------------------!
755! Description:
756! ------------
757!> @todo Missing subroutine description.
758!------------------------------------------------------------------------------!
759  SUBROUTINE interpolate_to_ocean( tag )
760
761#if defined( __parallel )
762
763    USE arrays_3d,                                                             &
764        ONLY:  total_2d_a, total_2d_o
765
766    USE indices,                                                               &
767        ONLY:  nbgp, nx, nx_a, nx_o, ny, ny_a, ny_o
768
769    USE kinds
770
771    USE pegrid
772
773    IMPLICIT NONE
774
775    INTEGER(iwp)             ::  dnx !<
776    INTEGER(iwp)             ::  dny !<
777    INTEGER(iwp)             ::  i   !<
778    INTEGER(iwp)             ::  ii  !<
779    INTEGER(iwp)             ::  j   !<
780    INTEGER(iwp)             ::  jj  !<
781    INTEGER(iwp), intent(in) ::  tag !<
782
783    REAL(wp)                 ::  fl  !<
784    REAL(wp)                 ::  fr  !<
785    REAL(wp)                 ::  myl !<
786    REAL(wp)                 ::  myr !<
787
788    CALL MPI_BARRIER( comm2d, ierr )
789
790    IF ( myid == 0 )  THEN   
791
792!
793!--    Number of gridpoints of the fine grid within one mesh of the coarse grid
794       dnx = ( nx_o + 1 ) / ( nx_a + 1 ) 
795       dny = ( ny_o + 1 ) / ( ny_a + 1 ) 
796
797!
798!--    Cyclic boundary conditions for atmosphere grid
799       total_2d_a(-nbgp:-1,:) = total_2d_a(ny+1-nbgp:ny,:)
800       total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx+1-nbgp:nx)
801       
802       total_2d_a(ny+1:ny+nbgp,:) = total_2d_a(0:nbgp-1,:)
803       total_2d_a(:,nx+1:nx+nbgp) = total_2d_a(:,0:nbgp-1)
804!
805!--    Bilinear Interpolation from atmosphere grid-layer to ocean grid-layer
806       DO  j = 0, ny
807          DO  i = 0, nx
808             myl = ( total_2d_a(j+1,i)   - total_2d_a(j,i)   ) / dny
809             myr = ( total_2d_a(j+1,i+1) - total_2d_a(j,i+1) ) / dny
810             DO  jj = 0, dny-1
811                fl = myl*jj + total_2d_a(j,i) 
812                fr = myr*jj + total_2d_a(j,i+1) 
813                DO  ii = 0, dnx-1
814                   total_2d_o(j*dny+jj,i*dnx+ii) = ( fr - fl ) / dnx * ii + fl
815                ENDDO
816             ENDDO
817          ENDDO
818       ENDDO
819!
820!--    Cyclic boundary conditions for ocean grid
821       total_2d_o(-nbgp:-1,:) = total_2d_o(ny_o+1-nbgp:ny_o,:)
822       total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx_o+1-nbgp:nx_o)
823
824       total_2d_o(ny_o+1:ny_o+nbgp,:) = total_2d_o(0:nbgp-1,:)
825       total_2d_o(:,nx_o+1:nx_o+nbgp) = total_2d_o(:,0:nbgp-1)
826
827       CALL MPI_SEND( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
828                      target_id, tag, comm_inter, ierr )
829
830    ENDIF
831
832    CALL MPI_BARRIER( comm2d, ierr ) 
833
834#endif
835
836  END SUBROUTINE interpolate_to_ocean
Note: See TracBrowser for help on using the repository browser.