Changeset 2232 for palm/trunk/SOURCE/surface_coupler.f90

Timestamp:

May 30, 2017 5:47:52 PM (7 years ago)

Author:

suehring

Message:

Adjustments according new topography and surface-modelling concept implemented

File:

• palm/trunk/SOURCE/surface_coupler.f90 (modified) (19 diffs)

palm/trunk/SOURCE/surface_coupler.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Adjust to new surface structure. Transfer 1D surface fluxes onto 2D grid 
+! (and back). 
 ! 
 ! Former revisions:
@@ -83 +84 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  pt, shf, qsws, qswst_remote, rho_ocean, sa, saswst, total_2d_a,       &
-               total_2d_o, tswst, u, usws, uswst, v, vsws, vswst
+        ONLY:  pt, rho_ocean, sa, total_2d_a, total_2d_o, u, v
 
     USE cloud_parameters,                                                      &
@@ -91 +91 @@
     USE control_parameters,                                                    &
         ONLY:  coupling_mode, coupling_mode_remote, coupling_topology,         &
-               humidity, humidity_remote, message_string, terminate_coupled,   &
-               terminate_coupled_remote, time_since_reference_point
+               humidity, humidity_remote, land_surface, message_string,        &
+               terminate_coupled, terminate_coupled_remote,                    &
+               time_since_reference_point, urban_surface
 
     USE cpulog,                                                                &
@@ -105 +106 @@
     USE pegrid
 
+    USE surface_mod,                                                           &
+        ONLY :  surf_def_h, surf_lsm_h, surf_type, surf_usm_h
+
     IMPLICIT NONE
 
+    INTEGER(iwp) ::  i                                    !< index variable x-direction 
+    INTEGER(iwp) ::  j                                    !< index variable y-direction 
+    INTEGER(iwp) ::  m                                    !< running index for surface elements
+
+    REAL(wp)    ::  cpw = 4218.0_wp                       !< heat capacity of water at constant pressure
     REAL(wp)    ::  time_since_reference_point_rem        !< 
     REAL(wp)    ::  total_2d(-nbgp:ny+nbgp,-nbgp:nx+nbgp) !< 
 
-    REAL(wp)    ::  cpw = 4218.0_wp !< heat capacity of water at constant pressure
+    REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  surface_flux !< dummy array for surface fluxes on 2D grid
+
 
 #if defined( __parallel )
@@ -164 +174 @@
 !
 !-- Exchange the current simulated time between the models,
-!-- currently just for total_2ding
+!-- currently just for total_2d
     IF ( coupling_topology == 0 ) THEN
    
@@ -197 +207 @@
 
 !
-!--       Send heat flux at bottom surface to the ocean
-          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
+!--       Send heat flux at bottom surface to the ocean. First, transfer from 
+!--       1D surface type to 2D grid.
+          CALL transfer_1D_to_2D_equal( surf_def_h(0)%shf, surf_lsm_h%shf,     &
+                                        surf_usm_h%shf )
+          CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, &
+                         12, comm_inter, ierr )
+!
+!--       Send humidity flux at bottom surface to the ocean. First, transfer
+!--       from 1D surface type to 2D grid.
+          CALL transfer_1D_to_2D_equal( surf_def_h(0)%qsws, surf_lsm_h%qsws,   &
+                                        surf_usm_h%qsws )
           IF ( humidity )  THEN
-             CALL MPI_SEND( qsws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 13, &
-                            comm_inter, ierr )
+             CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL,         &
+                            target_id, 13, comm_inter, ierr )
           ENDIF
 !
 !--       Receive temperature at the bottom surface from the ocean
-          CALL MPI_RECV( pt(0,nysg,nxlg), 1, type_xy, target_id, 14, &
+          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
-          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
-          CALL MPI_SEND( vsws(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
-                         comm_inter, ierr )
+!--       Send the momentum flux (u) at bottom surface to the ocean. First, 
+!--       transfer from 1D surface type to 2D grid.
+          CALL transfer_1D_to_2D_equal( surf_def_h(0)%usws, surf_lsm_h%usws,   &
+                                        surf_usm_h%usws )
+          CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, &
+                         15, comm_inter, ierr )
+!
+!--       Send the momentum flux (v) at bottom surface to the ocean. First, 
+!--       transfer from 1D surface type to 2D grid.
+          CALL transfer_1D_to_2D_equal( surf_def_h(0)%vsws, surf_lsm_h%vsws,   &
+                                        surf_usm_h%vsws )
+          CALL MPI_SEND( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, &
+                         16, comm_inter, ierr )
 !
 !--       Receive u at the bottom surface from the ocean
-          CALL MPI_RECV( u(0,nysg,nxlg), 1, type_xy, target_id, 17, &
+          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
-          CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, target_id, 18, &
+          CALL MPI_RECV( v(0,nysg,nxlg), 1, type_xy, target_id, 18,            &
                          comm_inter, status, ierr )
 !
@@ -235 +257 @@
           total_2d_a = 0.0_wp
           total_2d   = 0.0_wp
-          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, &
+!
+!--       Transfer from 1D surface type to 2D grid.
+          CALL transfer_1D_to_2D_unequal( surf_def_h(0)%shf, surf_lsm_h%shf,   &
+                                          surf_usm_h%shf )
+
+          CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0,  &
                            comm2d, ierr )
           CALL interpolate_to_ocean( 12 )   
@@ -245 +270 @@
              total_2d_a = 0.0_wp
              total_2d   = 0.0_wp
-             total_2d(nys:nyn,nxl:nxr) = qsws(nys:nyn,nxl:nxr)
+!
+!--          Transfer from 1D surface type to 2D grid.
+             CALL transfer_1D_to_2D_unequal( surf_def_h(0)%qsws,              &
+                                             surf_lsm_h%qsws,                 &
+                                             surf_usm_h%qsws )
 
              CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, &
@@ -254 +283 @@
 !--       Receive temperature at the bottom surface from the ocean
           IF ( myid == 0 )  THEN
-             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL, &
+             CALL MPI_RECV( total_2d_a(-nbgp,-nbgp), ngp_a, MPI_REAL,          &
                             target_id, 14, comm_inter, status, ierr )   
           ENDIF
@@ -265 +294 @@
           total_2d_a = 0.0_wp 
           total_2d   = 0.0_wp
-          total_2d(nys:nyn,nxl:nxr) = usws(nys:nyn,nxl:nxr)
+!
+!--       Transfer from 1D surface type to 2D grid.
+          CALL transfer_1D_to_2D_unequal( surf_def_h(0)%usws, surf_lsm_h%usws, &
+                                          surf_usm_h%usws )
           CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
                            comm2d, ierr )
@@ -273 +305 @@
           total_2d_a = 0.0_wp
           total_2d   = 0.0_wp
-          total_2d(nys:nyn,nxl:nxr) = vsws(nys:nyn,nxl:nxr)
+!
+!--       Transfer from 1D surface type to 2D grid.
+          CALL transfer_1D_to_2D_unequal( surf_def_h(0)%usws, surf_lsm_h%usws, &
+                                          surf_usm_h%usws )
           CALL MPI_REDUCE( total_2d, total_2d_a, ngp_a, MPI_REAL, MPI_SUM, 0, &
                            comm2d, ierr )
@@ -308 +343 @@
 !
 !--       Receive heat flux at the sea surface (top) from the atmosphere
-          CALL MPI_RECV( tswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
+          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 12, &
                          comm_inter, status, ierr )
+          CALL transfer_2D_to_1D_equal( surf_def_h(2)%shf )
 !
 !--       Receive humidity flux from the atmosphere (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, &
+             CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, &
                             target_id, 13, comm_inter, status, ierr )
+             CALL transfer_2D_to_1D_equal( surf_def_h(2)%qsws )
           ENDIF
 !
@@ -323 +360 @@
 !
 !--       Receive momentum flux (u) at the sea surface (top) from the atmosphere
-          CALL MPI_RECV( uswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
+          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 15, &
                          comm_inter, status, ierr )
+          CALL transfer_2D_to_1D_equal( surf_def_h(2)%usws )
 !
 !--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
-          CALL MPI_RECV( vswst(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
+          CALL MPI_RECV( surface_flux(nysg,nxlg), ngp_xy, MPI_REAL, target_id, 16, &
                          comm_inter, status, ierr )
+          CALL transfer_2D_to_1D_equal( surf_def_h(2)%vsws )
 !
 !--       Send u to the atmosphere
@@ -350 +389 @@
           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)
+          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%shf )
 !
 !--       Receive humidity flux at the sea surface (top) from the atmosphere
@@ -361 +400 @@
              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)
+             CALL transfer_2D_to_1D_unequal( surf_def_h(2)%qsws )
           ENDIF
 !
@@ -381 +420 @@
           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)
+          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%usws )
 !
 !--       Receive momentum flux (v) at the sea surface (top) from the atmosphere
@@ -391 +430 @@
           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)
+          CALL transfer_2D_to_1D_unequal( surf_def_h(2)%vsws )
 !
 !--       Send u to atmosphere
@@ -415 +454 @@
        IF ( humidity_remote )  THEN
 !
-!--       Here tswst is still the sum of atmospheric bottom heat fluxes,
+!--       Here top heat flux is still the sum of atmospheric bottom heat fluxes,
 !--       * latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
 !--       /(rho_atm(=1.0)*c_p)
-          tswst = tswst + qswst_remote * l_v / cp
-!
-!--        ...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_wp * sa(nzt,:,:) * 0.001_wp * qswst_remote /  &
-                    ( rho_ocean(nzt,:,:) * ( 1.0_wp - sa(nzt,:,:) * 0.001_wp ) )
+          DO  m = 1, surf_def_h(2)%ns
+             i = surf_def_h(2)%i(m)
+             j = surf_def_h(2)%j(m)
+             
+             surf_def_h(2)%shf(m) = surf_def_h(2)%shf(m) +                     &
+                                    surf_def_h(2)%qsws(m) * l_v / cp
+!
+!--          ...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 ) )
+             surf_def_h(2)%sasws(m) = -1.0_wp * sa(nzt,j,i) * 0.001_wp *       &
+                                      surf_def_h(2)%qsws(m) /                  &
+                                    ( rho_ocean(nzt,j,i) *                     &
+                                      ( 1.0_wp - sa(nzt,j,i) * 0.001_wp )      &
+                                    )
+          ENDDO
        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_ocean(nzt,:,:) * cp / cpw
-
-!
-!--    Adjust the momentum fluxes with respect to ocean density
-       uswst = uswst / rho_ocean(nzt,:,:)
-       vswst = vswst / rho_ocean(nzt,:,:)
+!--    (constants are the specific heat capacities for air and water), as well 
+!--    as momentum fluxes 
+       DO  m = 1, surf_def_h(2)%ns
+          i = surf_def_h(2)%i(m)
+          j = surf_def_h(2)%j(m)
+          surf_def_h(2)%shf(m) = surf_def_h(2)%shf(m) / rho_ocean(nzt,j,i) *   &
+                                 cp / cpw
+
+          surf_def_h(2)%usws(m) = surf_def_h(2)%usws(m) / rho_ocean(nzt,j,i)
+          surf_def_h(2)%vsws(m) = surf_def_h(2)%vsws(m) / rho_ocean(nzt,j,i)
+       ENDDO
 
     ENDIF
@@ -447 +498 @@
 
 #endif
+
+     CONTAINS 
+
+!       Description:
+!------------------------------------------------------------------------------!
+!>      Data transfer from 1D surface-data type to 2D dummy array for equal 
+!>      grids in atmosphere and ocean.
+!------------------------------------------------------------------------------!
+        SUBROUTINE transfer_1D_to_2D_equal( def_1d, lsm_1d, usm_1d )
+
+           IMPLICIT NONE
+
+            INTEGER(iwp) ::  i   !< running index x
+            INTEGER(iwp) ::  j   !< running index y
+            INTEGER(iwp) ::  m   !< running index surface type
+
+            REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  def_1d !< 1D surface flux, default surfaces
+            REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  lsm_1d !< 1D surface flux, natural surfaces
+            REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  usm_1d !< 1D surface flux, urban surfaces
+!
+!--         Transfer surface flux at default surfaces to 2D grid 
+            DO  m = 1, surf_def_h(0)%ns
+               i = surf_def_h(0)%i(m)
+               j = surf_def_h(0)%j(m)
+               surface_flux(j,i) = def_1d(m)
+            ENDDO
+!
+!--         Transfer surface flux at natural surfaces to 2D grid 
+            IF ( land_surface )  THEN 
+               DO  m = 1, SIZE(lsm_1d)
+                  i = surf_lsm_h%i(m)
+                  j = surf_lsm_h%j(m)
+                  surface_flux(j,i) = lsm_1d(m)
+               ENDDO
+            ENDIF
+!
+!--         Transfer surface flux at natural surfaces to 2D grid 
+            IF ( urban_surface )  THEN 
+               DO  m = 1, SIZE(usm_1d)
+                  i = surf_usm_h%i(m)
+                  j = surf_usm_h%j(m)
+                  surface_flux(j,i) = usm_1d(m)
+               ENDDO
+            ENDIF
+
+        END SUBROUTINE transfer_1D_to_2D_equal
+
+!       Description:
+!------------------------------------------------------------------------------!
+!>      Data transfer from 2D array for equal grids onto 1D surface-data type
+!>      array.
+!------------------------------------------------------------------------------!
+        SUBROUTINE transfer_2D_to_1D_equal( def_1d )
+
+           IMPLICIT NONE
+
+            INTEGER(iwp) ::  i   !< running index x
+            INTEGER(iwp) ::  j   !< running index y
+            INTEGER(iwp) ::  m   !< running index surface type
+
+            REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  def_1d !< 1D surface flux, default surfaces
+!
+!--         Transfer surface flux to 1D surface type, only for default surfaces 
+            DO  m = 1, surf_def_h(2)%ns
+               i = surf_def_h(2)%i(m)
+               j = surf_def_h(2)%j(m)
+               def_1d(m) = surface_flux(j,i)
+            ENDDO
+
+        END SUBROUTINE transfer_2D_to_1D_equal
+
+!       Description:
+!------------------------------------------------------------------------------!
+!>      Data transfer from 1D surface-data type to 2D dummy array from unequal 
+!>      grids in atmosphere and ocean.
+!------------------------------------------------------------------------------!
+        SUBROUTINE transfer_1D_to_2D_unequal( def_1d, lsm_1d, usm_1d )
+
+           IMPLICIT NONE
+
+            INTEGER(iwp) ::  i   !< running index x
+            INTEGER(iwp) ::  j   !< running index y
+            INTEGER(iwp) ::  m   !< running index surface type
+
+            REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  def_1d !< 1D surface flux, default surfaces
+            REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  lsm_1d !< 1D surface flux, natural surfaces
+            REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  usm_1d !< 1D surface flux, urban surfaces
+!
+!--         Transfer surface flux at default surfaces to 2D grid. Transfer no
+!--         ghost-grid points since total_2d is a global array.
+            DO  m = 1, SIZE(def_1d)
+               i = surf_def_h(0)%i(m)
+               j = surf_def_h(0)%j(m)
+
+               IF ( i >= nxl  .AND.  i <= nxr  .AND.                           &
+                    j >= nys  .AND.  j <= nyn )  THEN
+                  total_2d(j,i) = def_1d(m)
+               ENDIF
+            ENDDO
+!
+!--         Transfer surface flux at natural surfaces to 2D grid 
+            IF ( land_surface )  THEN 
+               DO  m = 1, SIZE(lsm_1d)
+                  i = surf_lsm_h%i(m)
+                  j = surf_lsm_h%j(m)
+
+                  IF ( i >= nxl  .AND.  i <= nxr  .AND.                        &
+                       j >= nys  .AND.  j <= nyn )  THEN
+                     total_2d(j,i) = lsm_1d(m)
+                  ENDIF
+               ENDDO
+            ENDIF
+!
+!--         Transfer surface flux at natural surfaces to 2D grid 
+            IF ( urban_surface )  THEN 
+               DO  m = 1, SIZE(usm_1d)
+                  i = surf_usm_h%i(m)
+                  j = surf_usm_h%j(m)
+
+                  IF ( i >= nxl  .AND.  i <= nxr  .AND.                        &
+                       j >= nys  .AND.  j <= nyn )  THEN
+                     total_2d(j,i) = usm_1d(m)
+                  ENDIF
+               ENDDO
+            ENDIF
+
+        END SUBROUTINE transfer_1D_to_2D_unequal
+
+!       Description:
+!------------------------------------------------------------------------------!
+!>      Data transfer from 2D dummy array from unequal grids to 1D surface-data 
+!>      type.
+!------------------------------------------------------------------------------!
+        SUBROUTINE transfer_2D_to_1D_unequal( def_1d )
+
+           IMPLICIT NONE
+
+            INTEGER(iwp) ::  i   !< running index x
+            INTEGER(iwp) ::  j   !< running index y
+            INTEGER(iwp) ::  m   !< running index surface type
+
+            REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  def_1d !< 1D surface flux, default surfaces
+!
+!--         Transfer 2D surface flux to default surfaces data type. Transfer no
+!--         ghost-grid points since total_2d is a global array.
+            DO  m = 1, SIZE(def_1d)
+               i = surf_def_h(2)%i(m)
+               j = surf_def_h(2)%j(m)
+
+               IF ( i >= nxl  .AND.  i <= nxr  .AND.                           &
+                    j >= nys  .AND.  j <= nyn )  THEN
+                  def_1d(m) = total_2d_o(j,i)
+               ENDIF
+            ENDDO
+
+
+        END SUBROUTINE transfer_2D_to_1D_unequal
 
   END SUBROUTINE surface_coupler