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

Last change on this file since 4785 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

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