Commit 96b20e92 authored by bjb's avatar bjb
Browse files

Change An distribution to two-field

parent 14f7834d
...@@ -116,6 +116,9 @@ ...@@ -116,6 +116,9 @@
allocate(An(E2DFIELD),stat=rc) allocate(An(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (An)' if (rc /= 0) stop 'init_2d: Error allocating memory (An)'
allocate(AnX(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (AnX)'
allocate(surfdiv(E2DFIELD),stat=rc) allocate(surfdiv(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_2d: Error allocating memory (surfdiv)' if (rc /= 0) stop 'init_2d: Error allocating memory (surfdiv)'
......
...@@ -21,7 +21,7 @@ ...@@ -21,7 +21,7 @@
REALTYPE,dimension(:,:),allocatable :: Slru,Slrv REALTYPE,dimension(:,:),allocatable :: Slru,Slrv
REALTYPE,dimension(:,:),allocatable :: zub,zvb REALTYPE,dimension(:,:),allocatable :: zub,zvb
REALTYPE,dimension(:,:),allocatable :: zub0,zvb0 REALTYPE,dimension(:,:),allocatable :: zub0,zvb0
REALTYPE,dimension(:,:),allocatable :: An REALTYPE,dimension(:,:),allocatable :: An,AnX
REALTYPE,dimension(:,:),allocatable :: surfdiv REALTYPE,dimension(:,:),allocatable :: surfdiv
REALTYPE,dimension(:,:),allocatable :: fwf,fwf_int REALTYPE,dimension(:,:),allocatable :: fwf,fwf_int
......
!$Id: m2d.F90,v 1.30 2009-05-13 09:46:10 bjb Exp $ !$Id: m2d.F90,v 1.31 2009-05-15 07:00:16 bjb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
use domain, only: ilg,ihg,jlg,jhg use domain, only: ilg,ihg,jlg,jhg
use domain, only: ill,ihl,jll,jhl use domain, only: ill,ihl,jll,jhl
use domain, only: openbdy,z0_method,z0_const,z0 use domain, only: openbdy,z0_method,z0_const,z0
use domain, only: az,ax
use halo_zones, only : update_2d_halo,wait_halo use halo_zones, only : update_2d_halo,wait_halo
use halo_zones, only : U_TAG,V_TAG,H_TAG use halo_zones, only : U_TAG,V_TAG,H_TAG
use variables_2d use variables_2d
...@@ -63,6 +64,8 @@ ...@@ -63,6 +64,8 @@
! Original author(s): Karsten Bolding & Hans Burchard ! Original author(s): Karsten Bolding & Hans Burchard
! !
! !LOCAL VARIABLES: ! !LOCAL VARIABLES:
integer :: num_neighbors
REALTYPE :: An_sum
! !
!EOP !EOP
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
...@@ -142,14 +145,50 @@ ...@@ -142,14 +145,50 @@
call getm_error("init_2d()", & call getm_error("init_2d()", &
"Constant horizontal numerical diffusion <0"); "Constant horizontal numerical diffusion <0");
else else
An = An_const An = An_const
AnX = An_const
end if end if
case(2) case(2)
LEVEL2 'An_method=2 -> Using space varying horizontal numerical diffusion' LEVEL2 'An_method=2 -> Using space varying horizontal numerical diffusion'
LEVEL2 '.. will read An from An_file ',trim(An_file) LEVEL2 '.. will read An from An_file ',trim(An_file)
! Initialize and then read field:
An = _ZERO_
call get_2d_field(trim(An_file),"An",ilg,ihg,jlg,jhg,An(ill:ihl,jll:jhl)) call get_2d_field(trim(An_file),"An",ilg,ihg,jlg,jhg,An(ill:ihl,jll:jhl))
call update_2d_halo(An,An,az,imin,jmin,imax,jmax,H_TAG) call update_2d_halo(An,An,az,imin,jmin,imax,jmax,H_TAG)
call wait_halo(H_TAG) call wait_halo(H_TAG)
! Compute AnX (An in X-points) based on An and the X- and T- masks
AnX = _ZERO_
do j=jmin,jmax+1
do i=imin,imax+1
if (ax(i,j) .ge. 1) then
num_neighbors = 0
An_sum = _ZERO_
! Each AnX should have up to 4 T-point neighbours.
if ( az(i-1,j-1) .ge. 1 ) then
An_sum = An_sum + An(i-1,j-1)
num_neighbors = num_neighbors +1
end if
if ( az(i-1,j) .ge. 1 ) then
An_sum = An_sum + An(i-1,j)
num_neighbors = num_neighbors +1
end if
if ( az(i,j-1) .ge. 1 ) then
An_sum = An_sum + An(i,j-1)
num_neighbors = num_neighbors +1
end if
if ( az(i,j) .ge. 1 ) then
An_sum = An_sum + An(i,j)
num_neighbors = num_neighbors +1
end if
! Take average of actual neighbours:
if (num_neighbors .gt. 0) then
AnX(i,j) = An_sum/num_neighbors
end if
end if
end do
end do
call update_2d_halo(AnX,AnX,ax,imin,jmin,imax,jmax,H_TAG)
call wait_halo(H_TAG)
case default case default
call getm_error("init_2d()", & call getm_error("init_2d()", &
"A non valid An method has been chosen"); "A non valid An method has been chosen");
...@@ -282,7 +321,7 @@ ...@@ -282,7 +321,7 @@
#ifndef UV_ADV_DIRECT #ifndef UV_ADV_DIRECT
call uv_advect() call uv_advect()
if (Am .gt. _ZERO_ .or. An_method .gt. 0) then if (Am .gt. _ZERO_ .or. An_method .gt. 0) then
call uv_diffusion(Am,An_method,An) ! Has to be called after uv_advect. call uv_diffusion(Am,An_method,An,AnX) ! Has to be called after uv_advect.
end if end if
call mirror_bdy_2d(UEx,U_TAG) call mirror_bdy_2d(UEx,U_TAG)
call mirror_bdy_2d(VEx,V_TAG) call mirror_bdy_2d(VEx,V_TAG)
......
...@@ -42,6 +42,7 @@ ...@@ -42,6 +42,7 @@
REALTYPE zub0(E2DFIELD) REALTYPE zub0(E2DFIELD)
REALTYPE zvb0(E2DFIELD) REALTYPE zvb0(E2DFIELD)
REALTYPE An(E2DFIELD) REALTYPE An(E2DFIELD)
REALTYPE AnX(E2DFIELD)
REALTYPE surfdiv(E2DFIELD) REALTYPE surfdiv(E2DFIELD)
REALTYPE fwf(E2DFIELD) REALTYPE fwf(E2DFIELD)
REALTYPE fwf_int(E2DFIELD) REALTYPE fwf_int(E2DFIELD)
......
!$Id: uv_diffusion.F90,v 1.10 2009-05-07 16:00:26 kb Exp $ !$Id: uv_diffusion.F90,v 1.11 2009-05-15 07:00:16 bjb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
! !ROUTINE: uv_diffusion - 2D diffusion of momentum \label{sec-uv-diffusion} ! !ROUTINE: uv_diffusion - 2D diffusion of momentum \label{sec-uv-diffusion}
! !
! !INTERFACE: ! !INTERFACE:
subroutine uv_diffusion(Am,An_method,An) subroutine uv_diffusion(Am,An_method,An,AnX)
! !
! !DESCRIPTION: ! !DESCRIPTION:
! !
...@@ -168,7 +168,7 @@ ...@@ -168,7 +168,7 @@
! !INPUT PARAMETERS: ! !INPUT PARAMETERS:
REALTYPE, intent(in) :: Am REALTYPE, intent(in) :: Am
integer, intent(in) :: An_method integer, intent(in) :: An_method
REALTYPE, intent(in) :: An(E2DFIELD) REALTYPE, intent(in) :: An(E2DFIELD),AnX(E2DFIELD)
! !
! !INPUT/OUTPUT PARAMETERS: ! !INPUT/OUTPUT PARAMETERS:
! !
...@@ -223,7 +223,7 @@ ...@@ -223,7 +223,7 @@
+(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX ) +(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX )
end if end if
if (An_method .gt. 0) then if (An_method .gt. 0) then
PP(i,j)=PP(i,j)+An(i,j)*(U(i,j+1)-U(i,j))*DXX/DYX PP(i,j)=PP(i,j)+AnX(i,j)*(U(i,j+1)-U(i,j))*DXX/DYX
end if end if
end if end if
end do end do
...@@ -248,7 +248,7 @@ ...@@ -248,7 +248,7 @@
+(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX ) +(V(i+1,j)/DV(i+1,j)-V(i,j)/DV(i,j))/DXX )
end if end if
if (An_method .gt. 0) then if (An_method .gt. 0) then
PP(i,j)=PP(i,j)+An(i,j)*(V(i+1,j)-V(i,j))*DYX/DXX PP(i,j)=PP(i,j)+AnX(i,j)*(V(i+1,j)-V(i,j))*DYX/DXX
end if end if
end if end if
end do end do
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment