Commit 9514766f authored by Knut's avatar Knut

Merge branch 'bdy_2d' into bdy_3d

Conflicts:
	src/3d/m3d.F90
	src/3d/uu_momentum_3d.F90
parents 2c049354 2726d010
......@@ -612,6 +612,7 @@
<fabm_adv_hor>4</fabm_adv_hor>
<fabm_adv_ver>4</fabm_adv_ver>
<fabm_AH>-1.0</fabm_AH>
<fabm_surface_flux_file></fabm_surface_flux_file>
</getm_fabm_nml>
</getm_fabm>
</scenario>
......@@ -923,6 +923,8 @@
</options>
</element>
<element name="fabm_AH" type="float" label="constant horizontal diffusivity of pelagic fabm fields" unit="m²/s"/>
<element name="fabm_surface_flux_file" type="string" label="name of file with FABM state variable surface fluxes">
</element>
</element>
</element>
</element>
......@@ -25,6 +25,7 @@
use domain, only: ilg,ihg,jlg,jhg
use domain, only: ill,ihl,jll,jhl
use domain, only: have_boundaries
!KB use get_field, only: get_2d_field
use advection, only: init_advection,print_adv_settings,NOADV
use halo_zones, only: update_2d_halo,wait_halo,H_TAG
use variables_2d
......@@ -82,13 +83,6 @@
REALTYPE,dimension(E2DFIELD),intent(out) :: phydiss
end subroutine physical_dissipation
! Temporary interface (should be read from module):
subroutine get_2d_field(fn,varname,il,ih,jl,jh,break_on_missing,f)
character(len=*),intent(in) :: fn,varname
integer, intent(in) :: il,ih,jl,jh
logical, intent(in) :: break_on_missing
REALTYPE, intent(out) :: f(:,:)
end subroutine get_2d_field
end interface
!
! !PUBLIC DATA MEMBERS:
......@@ -332,7 +326,7 @@
! of a hotstart file.
!
! !LOCAL VARIABLES:
integer :: i,j, ischange, rc
integer :: i,j, ischange
!EOP
!-------------------------------------------------------------------------
!BOC
......@@ -346,20 +340,7 @@
ufirst = ( mod(MinN,2) .eq. 0 )
if (do_numerical_analyses_2d) then
allocate(phydis_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (phydis_2d)'
phydis_2d = _ZERO_
allocate(numdis_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_2d)'
numdis_2d = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(numdis_2d_old(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_2d_old)'
numdis_2d_old = _ZERO_
#endif
end if
call postinit_variables_2d()
!
! It is possible that a user changes the land mask and reads an "old" hotstart file.
! In this case the "old" velocities will need to be zeroed out.
......
......@@ -159,6 +159,52 @@
end subroutine init_variables_2d
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: postinit_variables_2d - re-initialise some 2D stuff.
!
! !INTERFACE:
subroutine postinit_variables_2d()
IMPLICIT NONE
!
! !DESCRIPTION:
!
! !LOCAL VARIABLES:
integer :: rc
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'postinit_variables_2d() # ',Ncall
#endif
if (do_numerical_analyses_2d) then
allocate(phydis_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (phydis_2d)'
phydis_2d = _ZERO_
allocate(numdis_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_2d)'
numdis_2d = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(numdis_2d_old(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_2d_old)'
numdis_2d_old = _ZERO_
#endif
end if
#ifdef DEBUG
write(debug,*) 'Leaving postinit_variables_2d()'
write(debug,*)
#endif
return
end subroutine postinit_variables_2d
!EOC
!-----------------------------------------------------------------------
!BOP
!
......
......@@ -50,9 +50,9 @@
! onto the adapted grid afterwards.
!
! !USES:
use domain, only: ga,imin,imax,jmin,jmax,kmax,H,HU,HV,az,au,av
use domain, only: ga,imin,imax,jmin,jmax,kmax,H,az,au,av
use variables_3d, only: dt,kmin,kumin,kvmin,ho,hn,huo,hvo,hun,hvn
use variables_3d, only: Dn,Dun,Dvn,sseo,ssen,ssuo,ssvo
use variables_3d, only: Dn,Dun,Dvn,sseo,ssen
use variables_3d, only: kmin_pmz,kumin_pmz,kvmin_pmz
use variables_3d, only: preadapt
use vertical_coordinates,only: hcheck
......@@ -126,9 +126,7 @@
REALTYPE :: cNN,cSS,cdd,csum
REALTYPE :: cbg=6*_TENTH_
REALTYPE :: tfac_hor=3600*_ONE_ ! factor introducing a hor. adaption timescale = dt/tgrid_hor where tgrid_hor=tgrid/N
integer :: iip
integer,save :: count=0
namelist /adapt_coord/ faclag,facdif,fachor,faciso, &
depthmin,Ncrit, &
cNN,cSS,cdd,cbg,d_vel,d_dens, &
......
......@@ -56,8 +56,10 @@
#endif
if (first) then
if (.not. allocated(ga)) allocate(ga(0:kmax),stat=rc)
if (rc /= 0) stop 'coordinates: Error allocating (ga)'
if (.not. allocated(ga)) then
allocate(ga(0:kmax),stat=rc)
if (rc /= 0) stop 'coordinates: Error allocating (ga)'
end if
do k=0,kmax
ga(k) = k
end do
......
......@@ -16,6 +16,7 @@
use domain, only: ilg,ihg,jlg,jhg,ill,ihl,jll,jhl
use domain, only: az,latc,lonc
use domain,only: H
!KB use get_field, only: get_2d_field,get_3d_field
use variables_3d, only: uu,vv,ww,hun,hvn,ho,hn
use variables_3d,only: fabm_pel,fabm_ben,fabm_diag,fabm_diag_hz
use variables_3d, only: nuh,T,S,rho,a,g1,g2,taub
......@@ -23,27 +24,17 @@
use advection_3d, only: print_adv_settings_3d,do_advection_3d
use variables_2d, only: D,fwf_int
use meteo, only: swr,u10,v10,evap,precip,tcc
use time, only: yearday,secondsofday
use time, only: month,yearday,secondsofday,timestr
use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
! JORN_FABM
use gotm_fabm, only: init_gotm_fabm,set_env_gotm_fabm,do_gotm_fabm
use gotm_fabm, only: gotm_fabm_calc=>fabm_calc, model, cc_col=>cc, cc_diag_col=>cc_diag, cc_diag_hz_col=>cc_diag_hz, cc_transport
use fabm, only: type_horizontal_variable_id
use fabm, only: type_horizontal_variable_id, fabm_is_variable_used
use fabm_types,only: output_instantaneous, output_none
use fabm_standard_variables, only: standard_variables
IMPLICIT NONE
! Temporary interface (should be read from module):
interface
subroutine get_2d_field(fn,varname,il,ih,jl,jh,break_on_missing,f)
character(len=*),intent(in) :: fn,varname
integer, intent(in) :: il,ih,jl,jh
logical, intent(in) :: break_on_missing
REALTYPE, intent(out) :: f(:,:)
end subroutine get_2d_field
end interface
!
! !PUBLIC DATA MEMBERS:
public init_getm_fabm, postinit_getm_fabm, do_getm_fabm, model, output_none
......@@ -64,6 +55,21 @@
integer :: fabm_adv_ver=1
REALTYPE :: fabm_AH=-_ONE_
type (type_horizontal_variable_id) :: id_bottom_depth_below_geoid,id_bottom_depth
type type_input_variable
integer :: ncid = -1
integer :: varid = -1
class (type_input_variable), pointer :: next => null()
end type
type,extends(type_input_variable) :: type_horizontal_input_variable
REALTYPE, allocatable, dimension(:,:) :: data
type (type_horizontal_variable_id) :: id
end type
class (type_input_variable), pointer, save :: first_input_variable => null()
integer :: old_month=-1
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
......@@ -71,6 +77,23 @@
!EOP
!-----------------------------------------------------------------------
interface
subroutine inquire_file(fn,ncid,varids,varnames)
character(len=*), intent(in) :: fn
integer, intent(inout) :: ncid
integer, allocatable, intent(inout) :: varids(:)
character(len=50), allocatable, intent(out) :: varnames(:)
end subroutine inquire_file
!KB - only until a proper input_manager has been made
subroutine get_2d_field_ncdf_by_id(ncid,varid,il,ih,jl,jh,n,field)
integer, intent(in) :: ncid,varid
integer, intent(in) :: il,ih,jl,jh,n
REALTYPE, intent(out) :: field(:,:)
end subroutine get_2d_field_ncdf_by_id
end interface
contains
!-----------------------------------------------------------------------
......@@ -99,11 +122,16 @@
integer, parameter :: unit_fabm=63
integer :: rc
integer :: i,j,n
character(len=PATH_MAX) :: fabm_init_file
character(len=PATH_MAX) :: fabm_init_file, fabm_surface_flux_file=""
integer :: fabm_init_format, fabm_field_no
integer :: ncid
integer, allocatable :: varids(:)
character(len=50), allocatable :: varnames(:)
namelist /getm_fabm_nml/ fabm_init_method, &
fabm_init_file,fabm_init_format,fabm_field_no, &
fabm_surface_flux_file, &
fabm_adv_split,fabm_adv_hor,fabm_adv_ver,fabm_AH
!EOP
!-------------------------------------------------------------------------
......@@ -168,6 +196,26 @@
pa_nummix(n)%p3d => null()
end do
! Here we need to open the NetCDF file with FABM forcing data (if it exists)
! and loop over all its variables.
! For now 2D only, so make sure each NetCDF variable is 2D.
! For each variable, register_horizontal_input_variable should be called (see below).
! That looks up the FABM variable and also allocates the 2D field that will hold the input data.
LEVEL2 'FABM input and forcing ...'
if (len_trim(fabm_surface_flux_file) .ne. 0) then
LEVEL3 'reading surface fluxes from:'
LEVEL4 trim(fabm_surface_flux_file)
call inquire_file(fabm_surface_flux_file,ncid,varids,varnames)
do n=1,size(varids)
if ( varids(n) .ne. -1) then
! remeber surface_flux model in fabm.yaml
LEVEL4 'inquiring: ',trim(varnames(n))//'_flux'
call register_horizontal_input_variable(trim(varnames(n))//'_flux',ncid,varids(n))
end if
end do
end if
! Initialize biogeochemical state variables.
select case (fabm_init_method)
case(0)
......@@ -220,6 +268,51 @@
end subroutine init_getm_fabm
!EOC
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: register_horizontal_input_variable
!
! !INTERFACE:
subroutine register_horizontal_input_variable(name,ncid,varid)
!
! !DESCRIPTION:
! Registers FABM horizontal fluxes (surface)
!
! !USES:
IMPLICIT NONE
!
! !INPUT PARAMETERS:
character(len=*),intent(in) :: name
integer, intent(in) :: ncid,varid
!
! !REVISION HISTORY:
! See the log for the module
!
! !LOCAL VARIABLES
class (type_horizontal_input_variable), pointer :: variable
!EOP
!-------------------------------------------------------------------------
!BOC
! Create the input variable and set associated data (FABM id,
! NetCDF id, 2D data field).
allocate(variable)
variable%id = model%get_horizontal_variable_id(name)
if (.not.fabm_is_variable_used(variable%id)) then
LEVEL2 'Prescribed input variable '//trim(name)//' is not used by FABM.'
stop 'register_horizontal_input_variable: unrecognized variable name'
end if
variable%ncid = ncid
variable%varid = varid
allocate(variable%data(I2DFIELD))
variable%data = _ZERO_
! Prepend to the list of inout variables.
variable%next => first_input_variable
first_input_variable => variable
end subroutine register_horizontal_input_variable
!EOC
!-----------------------------------------------------------------------
!BOP
!
......@@ -298,22 +391,42 @@
REALTYPE :: bioshade(1:kmax)
REALTYPE :: wind_speed,I_0,taub_nonnorm,cloud
REALTYPE :: z(1:kmax)
class (type_input_variable), pointer :: current_input_variable
integer :: ncid,varid
logical :: some_var_ok=.false.
!EOP
!-----------------------------------------------------------------------
!BOC
call tic(TIM_GETM_FABM)
! First update all input fields
if (month .ne. old_month) then
old_month = month
current_input_variable => first_input_variable
do while (associated(current_input_variable))
select type (current_input_variable)
class is (type_horizontal_input_variable)
ncid = current_input_variable%ncid
varid = current_input_variable%varid
if (ncid .gt. 0 .and. varid .gt. 0) then
some_var_ok = .true.
call get_2d_field_ncdf_by_id(ncid,varid,ilg,ihg,jlg,jhg,month, &
current_input_variable%data(ill:ihl,jll:jhl))
end if
end select
current_input_variable => current_input_variable%next
end do
if (some_var_ok) then
LEVEL3 timestr,': reading FABM surface fluxes ...',month
end if
end if
do n=1,size(model%state_variables)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (az(i,j) .eq. 1) then
! Note (KK): This would be the correct dilution if hn was
! corrected in start_macro. This also requires,
! that ho=hn is done in coordinates!
! fabm_pel(i,j,kmax,n) = fabm_pel(i,j,kmax,n)*(_ONE_-fwf_int(i,j)/ho(i,j,kmax))
fabm_pel(i,j,kmax,n) = fabm_pel(i,j,kmax,n)* &
( ho(i,j,kmax) / (ho(i,j,kmax)+fwf_int(i,j)) )
fabm_pel(i,j,kmax,n) = fabm_pel(i,j,kmax,n)*(_ONE_-fwf_int(i,j)/ho(i,j,kmax))
end if
end do
end do
......@@ -379,6 +492,16 @@
call model%link_horizontal_data(id_bottom_depth_below_geoid,H(i,j))
call model%link_horizontal_data(id_bottom_depth,D(i,j))
! Transfer prescribed input variables
current_input_variable => first_input_variable
do while (associated(current_input_variable))
select type (current_input_variable)
class is (type_horizontal_input_variable)
call model%link_horizontal_data(current_input_variable%id,current_input_variable%data(i,j))
end select
current_input_variable => current_input_variable%next
end do
! Update biogeochemical variable values.
call do_gotm_fabm(kmax)
......
......@@ -345,67 +345,7 @@
ufirst = ( mod(int(ceiling((_ONE_*MinN)/M)),2) .eq. 1 )
if (do_numerical_analyses_3d) then
allocate(phydis_3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phydis_3d)'
phydis_3d = _ZERO_
allocate(phydis_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phydis_int)'
phydis_int = _ZERO_
allocate(numdis_3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_3d)'
numdis_3d = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(numdis_3d_old(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_3d_old)'
numdis_3d_old = _ZERO_
allocate(numdis_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_int)'
numdis_int = _ZERO_
#endif
if (update_temp) then
allocate(phymix_T(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_T)'
phymix_T = _ZERO_
allocate(phymix_T_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_T_int)'
phymix_T_int = _ZERO_
allocate(nummix_T(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_T)'
nummix_T = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(nummix_T_old(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_T_old)'
nummix_T_old = _ZERO_
allocate(nummix_T_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_T_int)'
nummix_T_int = _ZERO_
#endif
end if
if (update_salt) then
allocate(phymix_S(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S)'
phymix_S = _ZERO_
allocate(phymix_S_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phymix_S_int)'
phymix_S_int = _ZERO_
allocate(nummix_S(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_S)'
nummix_S = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(nummix_S_old(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_S_old)'
nummix_S_old = _ZERO_
allocate(nummix_S_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (nummix_S_int)'
nummix_S_int = _ZERO_
#endif
end if
end if
call postinit_variables_3d(update_temp,update_salt)
! Hotstart fix - see postinit_2d
if (hotstart) then
......
......@@ -18,6 +18,7 @@
use exceptions
use domain, only: imin,jmin,imax,jmax,kmax,ioff,joff
use domain, only: H,az
!KB use get_field, only: get_3d_field
use variables_2d, only: fwf_int
use variables_3d, only: S,hn,kmin
use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
......@@ -269,15 +270,7 @@
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (az(i,j) .eq. 1) then
! Note (KK): This would be the correct dilution if hn was
! corrected in start_macro. This also requires,
! that ho=hn is done in coordinates!
! S(i,j,kmax) = S(i,j,kmax)*(_ONE_-fwf_int(i,j)/ho(i,j,kmax))
! Developers note:
! The parentheses are set to minimize truncation errors for fwf_int=0
S(i,j,kmax) = S(i,j,kmax)* &
( ho(i,j,kmax) / (ho(i,j,kmax)+fwf_int(i,j)) )
S(i,j,kmax) = S(i,j,kmax)*(_ONE_-fwf_int(i,j)/ho(i,j,kmax))
end if
end do
end do
......
......@@ -24,10 +24,10 @@
! the water depths in all T-, U- and V-points to get the layer thicknesses.
!
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax,H,HU,HV
use domain, only: imin,imax,jmin,jmax,kmax,H
use domain, only: ga,ddu,ddl
use variables_3d, only: kmin,kumin,kvmin,ho,hn,huo,hun,hvo,hvn
use variables_3d, only: Dn,Dun,Dvn,sseo,ssuo,ssvo
use variables_3d, only: kmin,kumin,kvmin,ho,hn,hun,hvn
use variables_3d, only: Dn,Dun,Dvn,sseo
use vertical_coordinates,only: restart_with_ho,restart_with_hn
IMPLICIT NONE
!
......@@ -55,8 +55,10 @@
if (first) then
if (.not. allocated(ga)) allocate(ga(0:kmax),stat=rc)
if (rc /= 0) stop 'coordinates: Error allocating (ga)'
if (.not. allocated(ga)) then
allocate(ga(0:kmax),stat=rc)
if (rc /= 0) stop 'coordinates: Error allocating (ga)'
end if
allocate(dga(0:kmax),stat=rc)
if (rc /= 0) STOP 'coordinates: Error allocating (dga)'
ga(0) = -_ONE_
......
......@@ -57,14 +57,12 @@
call update_2d_halo(fwf_int,fwf_int,az,imin,jmin,imax,jmax,z_TAG)
call wait_halo(z_TAG)
! Note (KK): This would be the correct handling for fwf_int, but it
! requires ho=hn in coordinates. See also necessary changes
! in do_salinity and do_getm_fabm!
! Do not update ssen, because true sseo is needed!
! hn(:,:,kmax) = hn(:,:,kmax) + fwf_int
hn(:,:,kmax) = hn(:,:,kmax) + fwf_int
ssen = ssen + fwf_int
do j=jmin-HALO,jmax+HALO ! Defining 'old' and 'new' sea surface
do i=imin-HALO,imax+HALO ! elevation for macro time step
! Note (KK): this sseo already includes rivers and fwf
sseo(i,j)=ssen(i,j)
ssen(i,j)=z(i,j)
ssevel(i,j) = _HALF_ * ( sseo(i,j) + ssen(i,j) )
......
......@@ -18,6 +18,9 @@
! !USES:
use exceptions
use domain, only: imin,jmin,imax,kmax,jmax,H,az,dry_z
use domain, only: ill,ihl,jll,jhl
use domain, only: ilg,ihg,jlg,jhg
!KB use get_field, only: get_3d_field
use variables_3d, only: T,rad,hn,kmin,A,g1,g2
use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
IMPLICIT NONE
......@@ -39,6 +42,10 @@
REALTYPE :: temp_AH=-_ONE_
integer :: attenuation_method=0,jerlov=1
character(len=PATH_MAX) :: attenuation_file="attenuation.nc"
integer :: ncid=-1,A_id,g1_id,g2_id
integer, allocatable :: varids(:)
character(len=50), allocatable :: varnames(:)
integer :: old_month=-1
REALTYPE :: A_const=0.58,g1_const=0.35,g2_const=23.0
REALTYPE :: swr_bot_refl_frac=-_ONE_
REALTYPE :: swr_min_bot_frac=0.01
......@@ -51,6 +58,22 @@
!EOP
!-----------------------------------------------------------------------
interface
subroutine inquire_file(fn,ncid,varids,varnames)
character(len=*), intent(in) :: fn
integer, intent(inout) :: ncid
integer, allocatable, intent(inout) :: varids(:)
character(len=50), allocatable, intent(out) :: varnames(:)
end subroutine inquire_file
!KB - only until a proper input_manager has been made
subroutine get_2d_field_ncdf_by_id(ncid,varid,il,ih,jl,jh,n,field)
integer, intent(in) :: ncid,varid
integer, intent(in) :: il,ih,jl,jh,n
REALTYPE, intent(out) :: field(:,:)
end subroutine get_2d_field_ncdf_by_id
end interface
contains
!-----------------------------------------------------------------------
......@@ -144,12 +167,20 @@
case (1)
LEVEL3 'reading attenuation coefficients from:'
LEVEL4 trim(attenuation_file)
call inquire_file(attenuation_file,ncid,varids,varnames)
do n=1,size(varids)
if(trim(varnames(n)) .eq. "A") A_id = varids(n)
if(trim(varnames(n)) .eq. "g1") g1_id = varids(n)
if(trim(varnames(n)) .eq. "g2") g2_id = varids(n)
end do
#if 0
LEVEL1 'WARNING: reading routine not coded yet'
LEVEL1 'WARNING: setting to jerlov=1'
A_const=0.58;g1_const=0.35;g2_const=23.0
A=A_const
g1=g1_const
g2=g2_const
#endif