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

Last change on this file since 4422 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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