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

Last change on this file since 4327 was 4182, checked in by scharf, 5 years ago
  • corrected "Former revisions" section
  • minor formatting in "Former revisions" section
  • added "Author" section
  • Property svn:keywords set to Id
File size: 30.3 KB
RevLine 
[1682]1!> @file surface_coupler.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[3655]17! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[258]20! Current revisions:
[1092]21! ------------------
[1321]22!
[3049]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: surface_coupler.f90 4182 2019-08-22 15:20:23Z Giersch $
[4182]27! Corrected "Former revisions" section
28!
29! 3655 2019-01-07 16:51:22Z knoop
[3274]30! Modularization of all bulk cloud physics code components
[1321]31!
[4182]32! 109 2007-08-28 15:26:47Z letzel
33! Initial revision
34!
[102]35! Description:
36! ------------
[1682]37!> Data exchange at the interface between coupled models
[102]38!------------------------------------------------------------------------------!
[1682]39 SUBROUTINE surface_coupler
40 
[102]41
[1320]42    USE arrays_3d,                                                             &
[2232]43        ONLY:  pt, rho_ocean, sa, total_2d_a, total_2d_o, u, v
[1320]44
[3274]45    USE basic_constants_and_equations_mod,                                     &
46        ONLY:  c_p, l_v
[1427]47
[1320]48    USE control_parameters,                                                    &
49        ONLY:  coupling_mode, coupling_mode_remote, coupling_topology,         &
[2232]50               humidity, humidity_remote, land_surface, message_string,        &
51               terminate_coupled, terminate_coupled_remote,                    &
52               time_since_reference_point, urban_surface
[1320]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
[102]63    USE pegrid
64
[2232]65    USE surface_mod,                                                           &
66        ONLY :  surf_def_h, surf_lsm_h, surf_type, surf_usm_h
67
[102]68    IMPLICIT NONE
69
[2232]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
[1682]75    REAL(wp)    ::  time_since_reference_point_rem        !<
76    REAL(wp)    ::  total_2d(-nbgp:ny+nbgp,-nbgp:nx+nbgp) !<
[102]77
[2232]78    REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  surface_flux !< dummy array for surface fluxes on 2D grid
[1427]79
[2232]80
[206]81#if defined( __parallel )
[102]82
[667]83    CALL cpu_log( log_point(39), 'surface_coupler', 'start' )
[102]84
[667]85
86
[102]87!
[108]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.
[667]94
95    IF ( coupling_topology == 0 ) THEN
[709]96       CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, target_id, &
97                          0,                                                   &
98                          terminate_coupled_remote, 1, MPI_INTEGER, target_id, &
[667]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
[709]108       CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
109                       ierr )
[667]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
[108]116    IF ( terminate_coupled_remote > 0 )  THEN
[3045]117       WRITE( message_string, * ) 'remote model "',                            &
118                                  TRIM( coupling_mode_remote ),                &
119                                  '" terminated',                              &
[3046]120                                  '&with terminate_coupled_remote = ',         &
[3045]121                                  terminate_coupled_remote,                    &
[3046]122                                  '&local model  "', TRIM( coupling_mode ),    &
[3045]123                                  '" has',                                     &
[3046]124                                  '&terminate_coupled = ',                     &
[667]125                                   terminate_coupled
[258]126       CALL message( 'surface_coupler', 'PA0310', 1, 2, 0, 6, 0 )
[108]127       RETURN
128    ENDIF
[667]129 
[291]130
[108]131!
132!-- Exchange the current simulated time between the models,
[2232]133!-- currently just for total_2d
[709]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 )
[667]140    ELSE
[709]141
[667]142       IF ( myid == 0 ) THEN
[709]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,        &
[667]147                         target_id, 11, comm_inter, status, ierr )
[709]148
[667]149       ENDIF
[709]150
151       CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
152                       ierr )
153
[667]154    ENDIF
[102]155
156!
157!-- Exchange the interface data
158    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
[667]159   
160!
[709]161!--    Horizontal grid size and number of processors is equal in ocean and
162!--    atmosphere
163       IF ( coupling_topology == 0 )  THEN
[102]164
165!
[2232]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 )
[102]172!
[2232]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 )
[667]177          IF ( humidity )  THEN
[2232]178             CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL,         &
179                            target_id, 13, comm_inter, ierr )
[667]180          ENDIF
181!
[709]182!--       Receive temperature at the bottom surface from the ocean
[2232]183          CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, target_id, 14,           &
[709]184                         comm_inter, status, ierr )
[108]185!
[2232]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 )
[102]192!
[2232]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 )
[102]199!
[709]200!--       Receive u at the bottom surface from the ocean
[2232]201          CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, target_id, 17,            &
[709]202                         comm_inter, status, ierr )
[667]203!
[709]204!--       Receive v at the bottom surface from the ocean
[2232]205          CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, target_id, 18,            &
[709]206                         comm_inter, status, ierr )
[667]207!
208!--    Horizontal grid size or number of processors differs between
209!--    ocean and atmosphere
210       ELSE
211     
212!
[709]213!--       Send heat flux at bottom surface to the ocean
[1353]214          total_2d_a = 0.0_wp
215          total_2d   = 0.0_wp
[2232]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 )
[709]220
[2232]221          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0,  &
[709]222                           comm2d, ierr )
223          CALL interpolate_to_ocean( 12 )   
[667]224!
[709]225!--       Send humidity flux at bottom surface to the ocean
226          IF ( humidity )  THEN
[1353]227             total_2d_a = 0.0_wp
228             total_2d   = 0.0_wp
[2232]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 )
[709]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 )
[667]238          ENDIF
239!
[709]240!--       Receive temperature at the bottom surface from the ocean
241          IF ( myid == 0 )  THEN
[2232]242             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL,          &
[667]243                            target_id, 14, comm_inter, status, ierr )   
244          ENDIF
245          CALL MPI_BARRIER( comm2d, ierr )
[709]246          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
247                          ierr )
[667]248          pt(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
249!
[709]250!--       Send momentum flux (u) at bottom surface to the ocean
[1353]251          total_2d_a = 0.0_wp 
252          total_2d   = 0.0_wp
[2232]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 )
[709]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 )
[667]260!
[709]261!--       Send momentum flux (v) at bottom surface to the ocean
[1353]262          total_2d_a = 0.0_wp
263          total_2d   = 0.0_wp
[2232]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 )
[709]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 )
[667]271!
[709]272!--       Receive u at the bottom surface from the ocean
273          IF ( myid == 0 )  THEN
[667]274             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
[709]275                            target_id, 17, comm_inter, status, ierr )
[667]276          ENDIF
277          CALL MPI_BARRIER( comm2d, ierr )
[709]278          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
279                          ierr )
[667]280          u(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
281!
[709]282!--       Receive v at the bottom surface from the ocean
283          IF ( myid == 0 )  THEN
[667]284             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
[709]285                            target_id, 18, comm_inter, status, ierr )
[667]286          ENDIF
287          CALL MPI_BARRIER( comm2d, ierr )
[709]288          CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, 0, comm2d, &
289                          ierr )
[667]290          v(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg)
291
292       ENDIF
293
[102]294    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
295
296!
[667]297!--    Horizontal grid size and number of processors is equal
298!--    in ocean and atmosphere
299       IF ( coupling_topology == 0 ) THEN
300!
[709]301!--       Receive heat flux at the sea surface (top) from the atmosphere
[2232]302          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
[709]303                         comm_inter, status, ierr )
[2232]304          CALL transfer_2D_to_1D_equal( surf_def_h(2)%shf )
[102]305!
[709]306!--       Receive humidity flux from the atmosphere (bottom)
[667]307!--       and add it to the heat flux at the sea surface (top)...
308          IF ( humidity_remote )  THEN
[2232]309             CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, &
[667]310                            target_id, 13, comm_inter, status, ierr )
[2232]311             CALL transfer_2D_to_1D_equal( surf_def_h(2)%qsws )
[667]312          ENDIF
313!
314!--       Send sea surface temperature to the atmosphere model
[709]315          CALL MPI_SEND( pt(nzt,nysg,nxlg), 1, type_xy, target_id, 14, &
316                         comm_inter, ierr )
[667]317!
318!--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
[2232]319          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
[709]320                         comm_inter, status, ierr )
[2232]321          CALL transfer_2D_to_1D_equal( surf_def_h(2)%usws )
[667]322!
323!--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
[2232]324          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
[709]325                         comm_inter, status, ierr )
[2232]326          CALL transfer_2D_to_1D_equal( surf_def_h(2)%vsws )
[667]327!
[709]328!--       Send u to the atmosphere
329          CALL MPI_SEND( u(nzt,nysg,nxlg), 1, type_xy, target_id, 17, &
330                         comm_inter, ierr )
[667]331!
[709]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!
[667]336!--    Horizontal gridsize or number of processors differs between
337!--    ocean and atmosphere
338       ELSE
339!
[709]340!--       Receive heat flux at the sea surface (top) from the atmosphere
341          IF ( myid == 0 )  THEN
[667]342             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
[709]343                            target_id, 12, comm_inter, status, ierr )
[667]344          ENDIF
345          CALL MPI_BARRIER( comm2d, ierr )
[709]346          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
347                          ierr )
[2232]348          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%shf )
[667]349!
[709]350!--       Receive humidity flux at the sea surface (top) from the atmosphere
351          IF ( humidity_remote )  THEN
352             IF ( myid == 0 )  THEN
[667]353                CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
[709]354                               target_id, 13, comm_inter, status, ierr )
[667]355             ENDIF
356             CALL MPI_BARRIER( comm2d, ierr )
[709]357             CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, &
358                             comm2d, ierr)
[2232]359             CALL transfer_2D_to_1D_unequal( surf_def_h(2)%qsws )
[667]360          ENDIF
361!
362!--       Send surface temperature to atmosphere
[1353]363          total_2d_o = 0.0_wp
364          total_2d   = 0.0_wp
[667]365          total_2d(nys:nyn,nxl:nxr) = pt(nzt,nys:nyn,nxl:nxr)
366
[709]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 )
[667]370!
[709]371!--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
372          IF ( myid == 0 )  THEN
[667]373             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
[709]374                            target_id, 15, comm_inter, status, ierr )
[667]375          ENDIF
376          CALL MPI_BARRIER( comm2d, ierr )
377          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
[709]378                          0, comm2d, ierr )
[2232]379          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%usws )
[667]380!
[709]381!--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
382          IF ( myid == 0 )  THEN
[667]383             CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, &
[709]384                            target_id, 16, comm_inter, status, ierr )
[667]385          ENDIF
386          CALL MPI_BARRIER( comm2d, ierr )
[709]387          CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, 0, comm2d, &
388                          ierr )
[2232]389          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%vsws )
[667]390!
391!--       Send u to atmosphere
[1353]392          total_2d_o = 0.0_wp 
393          total_2d   = 0.0_wp
[667]394          total_2d(nys:nyn,nxl:nxr) = u(nzt,nys:nyn,nxl:nxr)
[709]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 )
[667]398!
399!--       Send v to atmosphere
[1353]400          total_2d_o = 0.0_wp
401          total_2d   = 0.0_wp
[667]402          total_2d(nys:nyn,nxl:nxr) = v(nzt,nys:nyn,nxl:nxr)
[709]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 )
[667]406       
407       ENDIF
408
409!
410!--    Conversions of fluxes received from atmosphere
411       IF ( humidity_remote )  THEN
[108]412!
[2232]413!--       Here top heat flux is still the sum of atmospheric bottom heat fluxes,
[709]414!--       * latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
415!--       /(rho_atm(=1.0)*c_p)
[2232]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) +                     &
[3274]421                                    surf_def_h(2)%qsws(m) * l_v / c_p
[709]422!
[2232]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
[108]432       ENDIF
433
434!
[102]435!--    Adjust the kinematic heat flux with respect to ocean density
[2232]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) *   &
[3274]442                                 c_p / cpw
[102]443
[2232]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
[102]447
[667]448    ENDIF
449
[709]450    IF ( coupling_topology == 1 )  THEN
[667]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
[2232]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
[667]615  END SUBROUTINE surface_coupler
616
617
618
[1682]619!------------------------------------------------------------------------------!
620! Description:
621! ------------
622!> @todo Missing subroutine description.
623!------------------------------------------------------------------------------!
[709]624  SUBROUTINE interpolate_to_atmos( tag )
[667]625
[880]626#if defined( __parallel )
627
[1320]628    USE arrays_3d,                                                             &
629        ONLY:  total_2d_a, total_2d_o
[667]630
[1320]631    USE indices,                                                               &
632        ONLY:  nbgp, nx, nx_a, nx_o, ny, ny_a, ny_o
633
634    USE kinds
635
[1324]636    USE pegrid
[1320]637
[667]638    IMPLICIT NONE
639
[1682]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   !<
[667]648
[1682]649    INTEGER(iwp), intent(in) ::  tag !<
[1320]650
[667]651    CALL MPI_BARRIER( comm2d, ierr )
652
[709]653    IF ( myid == 0 )  THEN
654!
655!--    Cyclic boundary conditions for the total 2D-grid
[667]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
[102]662!
[667]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) 
[102]666
667!
[709]668!--    Distance for interpolation around coarse grid points within the fine
669!--    grid (note: 2*dnx2 must not be equal with dnx)
[667]670       dnx2 = 2 * ( dnx / 2 )
671       dny2 = 2 * ( dny / 2 )
[102]672
[1353]673       total_2d_a = 0.0_wp
[102]674!
[667]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!
[709]688!--    Cyclic boundary conditions for atmosphere grid
[667]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
[709]696       CALL MPI_SEND( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, target_id, &
697                      tag, comm_inter, ierr )
[102]698
699    ENDIF
700
[667]701    CALL MPI_BARRIER( comm2d, ierr )
[102]702
[880]703#endif
704
[667]705  END SUBROUTINE interpolate_to_atmos
[102]706
[667]707
[1682]708!------------------------------------------------------------------------------!
709! Description:
710! ------------
711!> @todo Missing subroutine description.
712!------------------------------------------------------------------------------!
[709]713  SUBROUTINE interpolate_to_ocean( tag )
[667]714
[880]715#if defined( __parallel )
716
[1320]717    USE arrays_3d,                                                             &
718        ONLY:  total_2d_a, total_2d_o
[667]719
[1320]720    USE indices,                                                               &
721        ONLY:  nbgp, nx, nx_a, nx_o, ny, ny_a, ny_o
722
723    USE kinds
724
[1324]725    USE pegrid
[1320]726
[667]727    IMPLICIT NONE
728
[1682]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 !<
[667]736
[1682]737    REAL(wp)                 ::  fl  !<
738    REAL(wp)                 ::  fr  !<
739    REAL(wp)                 ::  myl !<
740    REAL(wp)                 ::  myr !<
[709]741
[667]742    CALL MPI_BARRIER( comm2d, ierr )
743
[709]744    IF ( myid == 0 )  THEN   
[667]745
746!
[709]747!--    Number of gridpoints of the fine grid within one mesh of the coarse grid
[667]748       dnx = ( nx_o + 1 ) / ( nx_a + 1 ) 
749       dny = ( ny_o + 1 ) / ( ny_a + 1 ) 
750
751!
[709]752!--    Cyclic boundary conditions for atmosphere grid
[667]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!
[709]759!--    Bilinear Interpolation from atmosphere grid-layer to ocean grid-layer
[667]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
[709]765                fl = myl*jj + total_2d_a(j,i) 
766                fr = myr*jj + total_2d_a(j,i+1) 
[667]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!
[709]774!--    Cyclic boundary conditions for ocean grid
[667]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
[880]788#endif
789
[667]790  END SUBROUTINE interpolate_to_ocean
Note: See TracBrowser for help on using the repository browser.