SUBROUTINE surface_coupler !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! additional case for nonequivalent processor and grid topopolgy in ocean and ! atmosphere added (coupling_topology = 1) ! ! ! Added exchange of u and v from Ocean to Atmosphere ! ! ! Former revisions: ! ------------------ ! $Id: surface_coupler.f90 667 2010-12-23 12:06:00Z suehring $ ! ! 291 2009-04-16 12:07:26Z raasch ! Coupling with independent precursor runs. ! Output of messages replaced by message handling routine. ! ! 206 2008-10-13 14:59:11Z raasch ! Implementation of a MPI-1 Coupling: replaced myid with target_id, ! deleted __mpi2 directives ! ! 109 2007-08-28 15:26:47Z letzel ! Initial revision ! ! Description: ! ------------ ! Data exchange at the interface between coupled models !------------------------------------------------------------------------------! USE arrays_3d USE control_parameters USE cpulog USE grid_variables USE indices USE interfaces USE pegrid IMPLICIT NONE INTEGER :: i, j, k REAL :: time_since_reference_point_rem REAL :: total_2d(-nbgp:ny+nbgp,-nbgp:nx+nbgp) #if defined( __parallel ) CALL cpu_log( log_point(39), 'surface_coupler', 'start' ) ! !-- In case of model termination initiated by the remote model !-- (terminate_coupled_remote > 0), initiate termination of the local model. !-- The rest of the coupler must then be skipped because it would cause an MPI !-- intercomminucation hang. !-- If necessary, the coupler will be called at the beginning of the next !-- restart run. IF ( coupling_topology == 0 ) THEN CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, target_id, & 0, & terminate_coupled_remote, 1, MPI_INTEGER, target_id, & 0, comm_inter, status, ierr ) ELSE IF ( myid == 0) THEN CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, & target_id, 0, & terminate_coupled_remote, 1, MPI_INTEGER, & target_id, 0, & comm_inter, status, ierr ) ENDIF CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr) ALLOCATE( total_2d_a(-nbgp:ny_a+nbgp,-nbgp:nx_a+nbgp), & total_2d_o(-nbgp:ny_o+nbgp,-nbgp:nx_o+nbgp) ) ENDIF IF ( terminate_coupled_remote > 0 ) THEN WRITE( message_string, * ) 'remote model "', & TRIM( coupling_mode_remote ), & '" terminated', & '&with terminate_coupled_remote = ', & terminate_coupled_remote, & '&local model "', TRIM( coupling_mode ), & '" has', & '&terminate_coupled = ', & terminate_coupled CALL message( 'surface_coupler', 'PA0310', 1, 2, 0, 6, 0 ) RETURN ENDIF ! !-- Exchange the current simulated time between the models, !-- currently just for total_2ding IF ( coupling_topology == 0 ) THEN CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, & target_id, 11, comm_inter, ierr ) CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, & target_id, 11, comm_inter, status, ierr ) ELSE IF ( myid == 0 ) THEN CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, & target_id, 11, comm_inter, ierr ) CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL, & target_id, 11, comm_inter, status, ierr ) ENDIF CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, & 0, comm2d, ierr ) ENDIF WRITE ( 9, * ) 'simulated time: ', simulated_time WRITE ( 9, * ) 'time since start of coupling: ', & time_since_reference_point, ' remote: ', & time_since_reference_point_rem CALL local_flush( 9 ) ! !-- Exchange the interface data IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN ! !-- Horizontal grid size and number of processors is equal !-- in ocean and atmosphere IF ( coupling_topology == 0 ) THEN ! !-- Send heat flux at bottom surface to the ocean model CALL MPI_SEND( shf(nysg,nxlg), ngp_xy, MPI_REAL, & target_id, 12, comm_inter, ierr ) ! !-- Send humidity flux at bottom surface to the ocean model IF ( humidity ) THEN CALL MPI_SEND( qsws(nysg,nxlg), ngp_xy, MPI_REAL, & target_id, 13, comm_inter, ierr ) ENDIF ! !-- Receive temperature at the bottom surface from the ocean model WRITE ( 9, * ) '*** receive pt from ocean' CALL local_flush( 9 ) CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, & target_id, 14, comm_inter, status, ierr ) ! !-- Send the momentum flux (u) at bottom surface to the ocean model CALL MPI_SEND( usws(nysg,nxlg), ngp_xy, MPI_REAL, & target_id, 15, comm_inter, ierr ) ! !-- Send the momentum flux (v) at bottom surface to the ocean model CALL MPI_SEND( vsws(nysg,nxlg), ngp_xy, MPI_REAL, & target_id, 16, comm_inter, ierr ) ! !-- Receive u at the bottom surface from the ocean model CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, & target_id, 17, comm_inter, status, ierr ) ! !-- Receive v at the bottom surface from the ocean model CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, & target_id, 18, comm_inter, status, ierr ) ! !-- Horizontal grid size or number of processors differs between !-- ocean and atmosphere ELSE ! !-- Send heat flux at bottom surface to the ocean model total_2d_a = 0.0 total_2d = 0.0 total_2d(nys:nyn,nxl:nxr) = shf(nys:nyn,nxl:nxr) CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, & MPI_SUM, 0, comm2d, ierr ) CALL interpolate_to_ocean(12) ! !-- Send humidity flux at bottom surface to the ocean model IF ( humidity ) THEN total_2d_a = 0.0 total_2d = 0.0 total_2d(nys:nyn,nxl:nxr) = qsws(nys:nyn,nxl:nxr) CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, & MPI_SUM, 0, comm2d, ierr ) CALL interpolate_to_ocean(13) ENDIF ! !-- Receive temperature at the bottom surface from the ocean model IF ( myid == 0 ) THEN CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, & target_id, 14, comm_inter, status, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, & 0, comm2d, ierr ) pt(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg) ! !-- Send momentum flux (u) at bottom surface to the ocean model total_2d_a = 0.0 total_2d = 0.0 total_2d(nys:nyn,nxl:nxr) = usws(nys:nyn,nxl:nxr) CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, & MPI_SUM, 0, comm2d, ierr ) CALL interpolate_to_ocean(15) ! !-- Send momentum flux (v) at bottom surface to the ocean model total_2d_a = 0.0 total_2d = 0.0 total_2d(nys:nyn,nxl:nxr) = vsws(nys:nyn,nxl:nxr) CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, & MPI_SUM, 0, comm2d, ierr ) CALL interpolate_to_ocean(16) ! !-- Receive u at the bottom surface from the ocean model IF ( myid == 0 ) THEN CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, & target_id, 17, comm_inter, status, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, & 0, comm2d, ierr ) u(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg) ! !-- Receive v at the bottom surface from the ocean model IF ( myid == 0 ) THEN CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, & target_id, 18, comm_inter, status, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) CALL MPI_BCAST( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, & 0, comm2d, ierr ) v(0,nysg:nyng,nxlg:nxrg) = total_2d_a(nysg:nyng,nxlg:nxrg) ENDIF ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN ! !-- Horizontal grid size and number of processors is equal !-- in ocean and atmosphere IF ( coupling_topology == 0 ) THEN ! !-- Receive heat flux at the sea surface (top) from the atmosphere model CALL MPI_RECV( tswst(nysg,nxlg), ngp_xy, MPI_REAL, & target_id, 12, comm_inter, status, ierr ) ! !-- Receive humidity flux from the atmosphere model (bottom) !-- and add it to the heat flux at the sea surface (top)... IF ( humidity_remote ) THEN CALL MPI_RECV( qswst_remote(nysg,nxlg), ngp_xy, MPI_REAL, & target_id, 13, comm_inter, status, ierr ) ENDIF ! !-- Send sea surface temperature to the atmosphere model CALL MPI_SEND( pt(nzt,nysg,nxlg), 1, type_xy, & target_id, 14, comm_inter, ierr ) ! !-- Receive momentum flux (u) at the sea surface (top) from the atmosphere !-- model WRITE ( 9, * ) '*** receive uswst from atmosphere' CALL local_flush( 9 ) CALL MPI_RECV( uswst(nysg,nxlg), ngp_xy, MPI_REAL, & target_id, 15, comm_inter, status, ierr ) ! !-- Receive momentum flux (v) at the sea surface (top) from the atmosphere !-- model CALL MPI_RECV( vswst(nysg,nxlg), ngp_xy, MPI_REAL, & target_id, 16, comm_inter, status, ierr ) !-- Send u to the atmosphere model CALL MPI_SEND( u(nzt,nysg,nxlg), 1, type_xy, & target_id, 17, comm_inter, ierr ) ! !-- Send v to the atmosphere model CALL MPI_SEND( v(nzt,nysg,nxlg), 1, type_xy, & target_id, 18, comm_inter, ierr ) ! !-- Horizontal gridsize or number of processors differs between !-- ocean and atmosphere ELSE ! !-- Receive heat flux at the sea surface (top) from the atmosphere model IF ( myid == 0 ) THEN CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & target_id, 12, comm_inter, status, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & 0, comm2d, ierr) tswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg) ! !-- Receive humidity flux at the sea surface (top) from the !-- atmosphere model IF ( humidity_remote ) THEN IF ( myid == 0 ) THEN CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & target_id, 13, comm_inter, status, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & 0, comm2d, ierr) qswst_remote(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg) ENDIF ! !-- Send surface temperature to atmosphere total_2d_o = 0.0 total_2d = 0.0 total_2d(nys:nyn,nxl:nxr) = pt(nzt,nys:nyn,nxl:nxr) CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, & MPI_REAL, MPI_SUM, 0, comm2d, ierr) CALL interpolate_to_atmos(14) ! !-- Receive momentum flux (u) at the sea surface (top) from the !-- atmosphere model IF ( myid == 0 ) THEN CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & target_id, 15, comm_inter, status, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & 0, comm2d, ierr) uswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg) ! !-- Receive momentum flux (v) at the sea surface (top) from the !-- atmosphere model IF ( myid == 0 ) THEN CALL MPI_RECV( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & target_id, 16, comm_inter, status, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) CALL MPI_BCAST( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & 0, comm2d, ierr) vswst(nysg:nyng,nxlg:nxrg) = total_2d_o(nysg:nyng,nxlg:nxrg) ! !-- Send u to atmosphere total_2d_o = 0.0 total_2d = 0.0 total_2d(nys:nyn,nxl:nxr) = u(nzt,nys:nyn,nxl:nxr) CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, MPI_REAL, & MPI_SUM, 0, comm2d, ierr) CALL interpolate_to_atmos(17) ! !-- Send v to atmosphere total_2d_o = 0.0 total_2d = 0.0 total_2d(nys:nyn,nxl:nxr) = v(nzt,nys:nyn,nxl:nxr) CALL MPI_REDUCE(total_2d, total_2d_o, ngp_o, MPI_REAL, & MPI_SUM, 0, comm2d, ierr) CALL interpolate_to_atmos(18) ENDIF ! !-- Conversions of fluxes received from atmosphere IF ( humidity_remote ) THEN !here tswst is still the sum of atmospheric bottom heat fluxes tswst = tswst + qswst_remote * 2.2626108e6 / 1005.0 !*latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol !/(rho_atm(=1.0)*c_p) ! !-- ...and convert it to a salinity flux at the sea surface (top) !-- following Steinhorn (1991), JPO 21, pp. 1681-1683: !-- S'w' = -S * evaporation / ( rho_water * ( 1 - S ) ) saswst = -1.0 * sa(nzt,:,:) * qswst_remote / & ( rho(nzt,:,:) * ( 1.0 - sa(nzt,:,:) ) ) ENDIF ! !-- Adjust the kinematic heat flux with respect to ocean density !-- (constants are the specific heat capacities for air and water) !-- now tswst is the ocean top heat flux tswst = tswst / rho(nzt,:,:) * 1005.0 / 4218.0 ! !-- Adjust the momentum fluxes with respect to ocean density uswst = uswst / rho(nzt,:,:) vswst = vswst / rho(nzt,:,:) ENDIF IF ( coupling_topology == 1 ) THEN DEALLOCATE( total_2d_o, total_2d_a ) ENDIF CALL cpu_log( log_point(39), 'surface_coupler', 'stop' ) #endif END SUBROUTINE surface_coupler SUBROUTINE interpolate_to_atmos(tag) USE arrays_3d USE control_parameters USE grid_variables USE indices USE pegrid IMPLICIT NONE INTEGER :: dnx, dnx2, dny, dny2, i, ii, j, jj INTEGER, intent(in) :: tag CALL MPI_BARRIER( comm2d, ierr ) IF ( myid == 0 ) THEN ! !-- cyclic boundary conditions for the total 2D-grid total_2d_o(-nbgp:-1,:) = total_2d_o(ny+1-nbgp:ny,:) total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx+1-nbgp:nx) total_2d_o(ny+1:ny+nbgp,:) = total_2d_o(0:nbgp-1,:) total_2d_o(:,nx+1:nx+nbgp) = total_2d_o(:,0:nbgp-1) ! !-- Number of gridpoints of the fine grid within one mesh of the coarse grid dnx = (nx_o+1) / (nx_a+1) dny = (ny_o+1) / (ny_a+1) ! !-- Distance for interpolation around coarse grid points within the fine grid !-- (note: 2*dnx2 must not be equal with dnx) dnx2 = 2 * ( dnx / 2 ) dny2 = 2 * ( dny / 2 ) total_2d_a = 0.0 ! !-- Interpolation from ocean-grid-layer to atmosphere-grid-layer DO j = 0, ny_a DO i = 0, nx_a DO jj = 0, dny2 DO ii = 0, dnx2 total_2d_a(j,i) = total_2d_a(j,i) & + total_2d_o(j*dny+jj,i*dnx+ii) ENDDO ENDDO total_2d_a(j,i) = total_2d_a(j,i) / ( ( dnx2 + 1 ) * ( dny2 + 1 ) ) ENDDO ENDDO ! !-- cyclic boundary conditions for atmosphere grid total_2d_a(-nbgp:-1,:) = total_2d_a(ny_a+1-nbgp:ny_a,:) total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx_a+1-nbgp:nx_a) total_2d_a(ny_a+1:ny_a+nbgp,:) = total_2d_a(0:nbgp-1,:) total_2d_a(:,nx_a+1:nx_a+nbgp) = total_2d_a(:,0:nbgp-1) ! !-- Transfer of the atmosphere-grid-layer to the atmosphere CALL MPI_SEND( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, & target_id, tag, comm_inter, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) END SUBROUTINE interpolate_to_atmos SUBROUTINE interpolate_to_ocean(tag) USE arrays_3d USE control_parameters USE grid_variables USE indices USE pegrid IMPLICIT NONE REAL :: fl, fr, myl, myr INTEGER :: dnx, dny, i, ii, j, jj INTEGER, intent(in) :: tag CALL MPI_BARRIER( comm2d, ierr ) IF ( myid == 0 ) THEN ! ! Number of gridpoints of the fine grid within one mesh of the coarse grid dnx = ( nx_o + 1 ) / ( nx_a + 1 ) dny = ( ny_o + 1 ) / ( ny_a + 1 ) ! !-- cyclic boundary conditions for atmosphere grid total_2d_a(-nbgp:-1,:) = total_2d_a(ny+1-nbgp:ny,:) total_2d_a(:,-nbgp:-1) = total_2d_a(:,nx+1-nbgp:nx) total_2d_a(ny+1:ny+nbgp,:) = total_2d_a(0:nbgp-1,:) total_2d_a(:,nx+1:nx+nbgp) = total_2d_a(:,0:nbgp-1) ! !-- Bilinear Interpolation from atmosphere-grid-layer to ocean-grid-layer DO j = 0, ny DO i = 0, nx myl = ( total_2d_a(j+1,i) - total_2d_a(j,i) ) / dny myr = ( total_2d_a(j+1,i+1) - total_2d_a(j,i+1) ) / dny DO jj = 0, dny-1 fl = myl*jj + total_2d_a(j,i) fr = myr*jj + total_2d_a(j,i+1) DO ii = 0, dnx-1 total_2d_o(j*dny+jj,i*dnx+ii) = ( fr - fl ) / dnx * ii + fl ENDDO ENDDO ENDDO ENDDO ! !-- cyclic boundary conditions for ocean grid total_2d_o(-nbgp:-1,:) = total_2d_o(ny_o+1-nbgp:ny_o,:) total_2d_o(:,-nbgp:-1) = total_2d_o(:,nx_o+1-nbgp:nx_o) total_2d_o(ny_o+1:ny_o+nbgp,:) = total_2d_o(0:nbgp-1,:) total_2d_o(:,nx_o+1:nx_o+nbgp) = total_2d_o(:,0:nbgp-1) CALL MPI_SEND( total_2d_o(-nbgp,-nbgp), ngp_o, MPI_REAL, & target_id, tag, comm_inter, ierr ) ENDIF CALL MPI_BARRIER( comm2d, ierr ) END SUBROUTINE interpolate_to_ocean