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

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

last commit documented

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