Skip to content

GitLab

Commit ea8effdc authored by Knut's avatar Knut
Browse files

finished integration of flexible output

parent 311b280d
Hide whitespace changes
Inline Side-by-side
Showing with 138 additions and 136 deletions
src/2d/depth_update.F90
View file @ ea8effdc
......@@ -44,6 +44,7 @@
CALL tic(TIM_DPTHUPDATE)
! TODO/BJB: Why is this turned off?
! KK: ...because we need to have non-zero DU/DV at land-sea-interfaces
#undef USE_MASK
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,d1,d2i,x)
......@@ -53,7 +54,6 @@
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
! TODO/BJB: Is it enough to do this on az?
D(i,j) = z(i,j)+H(i,j)
end do
end do
......
src/2d/m2d.F90
View file @ ea8effdc
......@@ -179,7 +179,7 @@
select case (elev_method)
case(1)
LEVEL2 'setting initial surface elevation to ',real(elev_const)
z = elev_const
where ( az.gt.0 ) z = elev_const
case(2)
LEVEL2 'getting initial surface elevation from ',trim(elev_file)
call get_2d_field(trim(elev_file),"elev",ilg,ihg,jlg,jhg,.true.,z(ill:ihl,jll:jhl))
......@@ -190,7 +190,7 @@
stop 'init_2d(): invalid elev_method'
end select
where ( z .lt. -H+min_depth)
where ( az.gt.0 .and. z.lt.-H+min_depth)
z = -H+min_depth
end where
zo = z
......@@ -387,13 +387,14 @@
V = _ZERO_
Vinto = _ZERO_
end where
! This is probably not absolutely necessary:
where (az .eq. 0)
z = _ZERO_
zo = _ZERO_
end where
end if
! This is only needed for proper flexible output
where (az .eq. 0)
z = -9999._rk
zo = -9999._rk
end where
call depth_update()
end if
......@@ -472,7 +473,7 @@
call toc(TIM_INTEGR2D)
end if
if (have_boundaries) call update_2d_bdy(loop,bdy2d_ramp)
call sealevel()
call sealevel(loop)
call depth_update()
if(residual .gt. 0) then
......
src/2d/sealevel.F90
View file @ ea8effdc
......@@ -5,7 +5,7 @@
! !ROUTINE: sealevel - using the cont. eq. to get the sealevel.
!
! !INTERFACE:
subroutine sealevel
subroutine sealevel(loop)
!
! !DESCRIPTION:
!
......@@ -26,7 +26,7 @@
#else
use domain, only : dx,dy,ard1
#endif
use m2d, only: dtm
use m2d, only: dtm,sealevel_check
use variables_2d, only: z,zo,U,V,fwf
use getm_timers, only: tic, toc, TIM_SEALEVEL, TIM_SEALEVELH
use halo_zones, only : update_2d_halo,wait_halo,z_TAG
......@@ -37,6 +37,9 @@
#endif
!$ use omp_lib
IMPLICIT NONE
! !INPUT PARAMETERS:
integer, intent(in) :: loop
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
......@@ -137,7 +140,9 @@
end do !end do while(break_flag>0)
#endif
call sealevel_nan_check()
if (sealevel_check.ne.0 .and. mod(loop,abs(sealevel_check)).eq.0) then
call sealevel_nan_check()
end if
#ifdef SLICE_MODEL
do i=imin,imax
......@@ -217,7 +222,7 @@
!
! Fast return if we should not sweep:
if (sealevel_check.eq.0 .or. have_warned.gt.0) then
if ( have_warned .gt. 0 ) then
#ifdef DEBUG
write(debug,*) 'Leaving sealevel_nan_check() early'
write(debug,*)
......@@ -263,7 +268,7 @@
end if
!
! Main bulk of work:
if (can_check .gt. 0 .and. mod(Ncall,abs(sealevel_check)).eq.0) then
if (can_check .gt. 0 ) then
! Count number of NaNs encountered
num_nan = 0
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,idum)
......
src/2d/variables_2d.F90
View file @ ea8effdc
......@@ -22,6 +22,7 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
integer, parameter :: rk = kind(_ONE_)
#ifdef STATIC
#include "static_2d.h"
#else
......@@ -33,8 +34,6 @@
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
!
! !LOCAL VARIABLES:
integer :: rc
!EOP
!-----------------------------------------------------------------------
......@@ -92,6 +91,8 @@
! !INPUT PARAMETERS:
integer, intent(in) :: runtype
!
! !LOCAL VARIABLES:
integer :: rc
!EOP
!-------------------------------------------------------------------------
!BOC
......@@ -114,10 +115,11 @@
break_stat = 0
#endif
z = _ZERO_; zo =_ZERO_
z = -9999._rk ; zo =_ZERO_
zub=_ZERO_ ; zub0=_ZERO_
zvb=_ZERO_ ; zvb0=_ZERO_
D = _ZERO_;
U = _ZERO_; DU = _ZERO_; fU = _ZERO_; Uint = _ZERO_; UEx = _ZERO_
V = _ZERO_; DV = _ZERO_; fV = _ZERO_; Vint = _ZERO_; VEx = _ZERO_
......@@ -168,7 +170,6 @@
!
! !LOCAL VARIABLES:
logical :: used
integer,parameter :: rk = kind(_ONE_)
!EOP
!-----------------------------------------------------------------------
!BOC
......
src/3d/getm_fabm.F90
View file @ ea8effdc
......@@ -292,10 +292,24 @@ end interface
! Update halos with biogeochemical variable values (distribute initial values).
do n=1,size(model%state_variables)
fabm_pel(:,:,0,n) = model%state_variables(n)%missing_value
forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) &
fabm_pel(i,j,:,n) = model%state_variables(n)%missing_value
call update_3d_halo(fabm_pel(:,:,:,n),fabm_pel(:,:,:,n),az, &
imin,jmin,imax,jmax,kmax,D_TAG)
call wait_halo(D_TAG)
end do
do n=1,size(model%bottom_state_variables)
where ( az.eq.0 ) fabm_ben(:,:,n) = model%bottom_state_variables(n)%missing_value
end do
do n=1,size(model%diagnostic_variables)
fabm_diag(:,:,0,n) = model%diagnostic_variables(n)%missing_value
forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) &
fabm_diag(i,j,:,n) = model%diagnostic_variables(n)%missing_value
end do
do n=1,size(model%horizontal_diagnostic_variables)
where ( az.eq.0 ) fabm_diag_hz(:,:,n) = model%horizontal_diagnostic_variables(n)%missing_value
end do
#ifdef DEBUG
write(debug,*) 'Leaving init_getm_fabm_fields()'
......@@ -366,24 +380,6 @@ end interface
fill_value=model%horizontal_diagnostic_variables(i)%missing_value, data2d=fabm_diag_hz(_2D_W_,i), category='fabm'//model%horizontal_diagnostic_variables(i)%target%owner%get_path(), output_level=output_level, used=in_output)
if (in_output) model%horizontal_diagnostic_variables(i)%save = .true.
end do
#if 0
if (do_numerical_analyses_3d) then
do i=1,size(model%state_variables)
call fm%register('nummix_'//trim(model%state_variables(i)%name), &
'('//trim(model%state_variables(i)%units)//')**2/s', &
'numerical mixing of '//trim(model%state_variables(i)%long_name), &
dimensions=(/id_dim_z/), &
category='fabm'//model%state_variables(i)%target%owner%get_path(), &
output_level=output_level_debug)
call fm%register('phymix_'//trim(model%state_variables(i)%name), &
'('//trim(model%state_variables(i)%units)//')**2/s', &
'physical mixing of '//trim(model%state_variables(i)%long_name), &
dimensions=(/id_dim_z/), &
category='fabm'//model%state_variables(i)%target%owner%get_path(), &
output_level=output_level_debug)
end do
end if
#endif
return
end subroutine register_fabm_variables
......
src/3d/m3d.F90
View file @ ea8effdc
......@@ -86,6 +86,7 @@
!
! !INTERFACE:
subroutine init_3d(runtype,timestep,hotstart)
!
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -93,7 +94,6 @@
REALTYPE, intent(in) :: timestep
logical, intent(in) :: hotstart
!
!
! !DESCRIPTION:
! Here, the {\tt m3d} namelist is read from {\tt getm.inp}, and the
! initialisation of variables is called (see routine {\tt init\_variables}
......@@ -397,8 +397,8 @@
num(i,j,:) = 1.e-15
nuh(i,j,:) = 1.e-15
#ifndef NO_BAROCLINIC
S(i,j,:) = _ZERO_
T(i,j,:) = _ZERO_
S(i,j,:) = -9999._rk
T(i,j,:) = -9999._rk
#endif
end if
end do
......
src/3d/salinity.F90
View file @ ea8effdc
......@@ -20,7 +20,7 @@
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 variables_3d, only: rk,S,hn,kmin
use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
IMPLICIT NONE
!
......@@ -123,19 +123,20 @@
LEVEL4 'out-of-bound values result in warnings only'
end if
if (salt_method .ne. 0) then
call check_3d_fields(imin,jmin,imax,jmax,kmin,kmax,az, &
S,min_salt,max_salt,status)
if (status .gt. 0) then
if (salt_check .gt. 0) then
call getm_error("do_salinity()", &
call getm_error("init_salinity()", &
"out-of-bound values encountered")
end if
if (salt_check .lt. 0) then
LEVEL1 'do_salinity(): ',status, &
LEVEL1 'init_salinity(): ',status, &
' out-of-bound values encountered'
end if
end if
end if
end if
......@@ -200,6 +201,9 @@
stop 'init_salinity'
end select
S(:,:,0) = -9999._rk
forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) S(i,j,:) = -9999._rk
call update_3d_halo(S,S,az,imin,jmin,imax,jmax,kmax,D_TAG)
call wait_halo(D_TAG)
call mirror_bdy_3d(S,D_TAG)
......
src/3d/temperature.F90
View file @ ea8effdc
......@@ -21,7 +21,7 @@
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 variables_3d, only: rk,T,rad,hn,kmin,A,g1,g2
use halo_zones, only: update_3d_halo,wait_halo,D_TAG,H_TAG
IMPLICIT NONE
!
......@@ -202,18 +202,20 @@ end interface
LEVEL4 'out-of-bound values result in warnings only'
end if
if (temp_method .ne. 0) then
call check_3d_fields(imin,jmin,imax,jmax,kmin,kmax,az, &
T,min_temp,max_temp,status)
if (status .gt. 0) then
if (temp_check .gt. 0) then
call getm_error("do_temperature()", &
call getm_error("init_temperature()", &
"out-of-bound values encountered")
end if
if (temp_check .lt. 0) then
LEVEL1 'do_temperature(): ',status, &
LEVEL1 'init_temperature(): ',status, &
' out-of-bound values encountered'
end if
end if
end if
end if
#ifdef DEBUG
......@@ -277,6 +279,9 @@ end interface
stop 'init_temperature'
end select
T(:,:,0) = -9999._rk
forall(i=imin:imax,j=jmin:jmax, az(i,j).eq.0) T(i,j,:) = -9999._rk
call update_3d_halo(T,T,az,imin,jmin,imax,jmax,kmax,D_TAG)
call wait_halo(D_TAG)
call mirror_bdy_3d(T,D_TAG)
......
src/3d/variables_3d.F90
View file @ ea8effdc
......@@ -116,6 +116,7 @@
IMPLICIT NONE
!
! !PUBLIC DATA MEMBERS:
integer, parameter :: rk = kind(_ONE_)
REALTYPE :: dt,cnpar=0.9
REALTYPE :: avmback=_ZERO_,avhback=_ZERO_
logical :: do_numerical_analyses=.false.
......@@ -192,8 +193,11 @@
#include "dynamic_allocations_3d.h"
#endif
kmin = 1
hn = _ZERO_ ; hun = _ZERO_ ; hvn = _ZERO_
uu = _ZERO_ ; vv = _ZERO_ ; ww = _ZERO_
#ifdef _MOMENTUM_TERMS_
tdv_u = _ZERO_ ; adv_u = _ZERO_ ; vsd_u = _ZERO_ ; hsd_u = _ZERO_
cor_u = _ZERO_ ; epg_u = _ZERO_ ; ipg_u = _ZERO_
......@@ -215,8 +219,8 @@
#endif
! must be nonzero for gotm_fabm in case of calc_temp=F
g1 = -9999*_ONE_
g2 = -9999*_ONE_
g1 = -9999._rk
g2 = -9999._rk
#ifdef DEBUG
write(debug,*) 'Leaving init_variables_3d()'
......@@ -249,7 +253,7 @@
! Original author(s): Karsten Bolding & Jorn Bruggeman
!
! !LOCAL VARIABLES:
integer,parameter :: rk = kind(_ONE_)
!
!EOP
!-----------------------------------------------------------------------
!BOC
......
src/futils/getm_version.F90
View file @ ea8effdc
module getm_version
character(len=*),parameter :: git_commit_id = "2.5.0"
character(len=*),parameter :: git_branch_name = "master"
character(len=*),parameter :: git_branch_name = "output"
end module
src/futils/time.F90
View file @ ea8effdc
......@@ -27,7 +27,7 @@
IMPLICIT NONE
!
! !PUBLIC MEMBER FUNCTIONS:
public :: time_diff
public :: time_diff,CalDat,JulDay
! !PUBLIC DATA MEMBERS:
integer, target :: julianday,secondsofday,yearday
......
src/getm/initialise.F90
View file @ ea8effdc
......@@ -351,9 +351,12 @@
end if
if (.not. dryrun) then
call do_output(runtype,MinN-1,timestep)
if (save_initial) then
call output_manager_prepare_save(julianday, int(secondsofday), 0, int(MinN-1))
call do_output_processing()
if (save_initial) call output_manager_save(julianday,secondsofday,MinN-1)
call output_manager_save(julianday,secondsofday,MinN-1)
end if
call do_output(runtype,MinN-1,timestep)
end if
#ifdef DEBUG
......
src/getm/integration.F90
View file @ ea8effdc
......@@ -62,8 +62,8 @@
use suspended_matter, only: spm_calc,do_spm
#endif
use input, only: do_input
use output, only: do_output,meanout
use output_processing, only: do_output_processing
use output, only: do_output
#ifdef TEST_NESTING
use nesting, only: nesting_file
#endif
......@@ -102,10 +102,6 @@
LEVEL1 t_(1:2),':',t_(3:4),':',t_(5:10),' n=',n
end if
call tic(TIM_FLEX_OUTPUT)
call output_manager_prepare_save(julianday, int(secondsofday), 0, int(n))
call toc(TIM_FLEX_OUTPUT)
#ifndef NO_3D
do_3d = (runtype .ge. 2 .and. mod(n,M) .eq. 0)
#endif
......@@ -157,19 +153,18 @@
call nesting_file(WRITING)
end if
#endif
call update_time(n)
call tic(TIM_FLEX_OUTPUT)
call output_manager_prepare_save(julianday, int(secondsofday), 0, int(n))
call toc(TIM_FLEX_OUTPUT)
call tic(TIM_OUTPUT_PROC)
call do_output_processing()
call toc(TIM_OUTPUT_PROC)
call tic(TIM_FLEX_OUTPUT)
call output_manager_save(julianday,secondsofday,n)
call toc(TIM_FLEX_OUTPUT)
call update_time(n)
#ifndef NO_3D
if(meanout .ge. 0) then
call calc_mean_fields(n,meanout)
end if
#endif
call do_output(runtype,n,timestep)
#ifdef DIAGNOSE
call diagnose(n,MaxN,runtype)
......
src/getm/register_all_variables.F90
View file @ ea8effdc
......@@ -354,26 +354,11 @@
! Original author(s): Karsten Bolding & Jorn Bruggeman
!
! !LOCAL VARIABLES:
integer :: i
!EOP
!-----------------------------------------------------------------------
!BOC
LEVEL1 'finalize_register_all_variables()'
! category - fabm
#if 0
#ifdef _FABM_
if (fabm_calc) then
if (do_numerical_analyses_3d) then
do i=1,size(model%state_variables)
call fm%send_data('nummix_'//trim(model%state_variables(i)%name), nummix_fabm_pel(_3D_W_,i))
call fm%send_data('phymix_'//trim(model%state_variables(i)%name), phymix_fabm_pel(_3D_W_,i))
end do
end if
end if
#endif
#endif
return
end subroutine finalize_register_all_variables
!EOC
......
src/meteo/meteo.F90
View file @ ea8effdc
......@@ -158,6 +158,7 @@
LEVEL3 'ty = ',ty
LEVEL3 'swr = ',swr_const
LEVEL3 'shf = ',shf_const
calc_met = .false.
case (2)
if(on_grid) then
LEVEL2 'Meteorological fields are on the computational grid'
......
src/ncdf/init_2d_ncdf.F90
View file @ ea8effdc
......@@ -46,11 +46,6 @@
! initialize all time-independent, grid related variables
call init_grid_ncdf(ncid,init3d,x_dim,y_dim)
! allocate workspace
allocate(ws(E2DFIELD),stat=err)
if (err .ne. 0) call getm_error("init_2d_ncdf()", &
"error allocating ws")
! define unlimited dimension
err = nf90_def_dim(ncid,'time',NF90_UNLIMITED,time_dim)
if (err .NE. NF90_NOERR) go to 10
......
src/ncdf/init_3d_ncdf.F90
View file @ ea8effdc
......@@ -13,7 +13,6 @@
use netcdf
use exceptions
use ncdf_common
use ncdf_2d, only: ws2d => ws
use ncdf_3d
use domain, only: ioff,joff
use domain, only: imin,imax,jmin,jmax,kmax
......@@ -56,16 +55,6 @@
! initialize all time-independent, grid related variables
call init_grid_ncdf(ncid,init3d,x_dim,y_dim,z_dim)
! allocate workspace
if (.not. allocated(ws2d)) then
allocate(ws2d(E2DFIELD),stat=err)
if (err .ne. 0) call getm_error("init_3d_ncdf()", &
"error allocating ws2d")
end if
allocate(ws(I3DFIELD),stat=err)
if (err .ne. 0) call getm_error("init_3d_ncdf()", &
"error allocating ws")
! define unlimited dimension
err = nf90_def_dim(ncid,'time',NF90_UNLIMITED,time_dim)
if (err .NE. NF90_NOERR) go to 10
......@@ -283,41 +272,40 @@
#endif
end if
if (save_strho) then
if (calc_salt .and. save_s) then
fv = salt_missing
mv = salt_missing
vr(1) = 0.
vr(2) = 40.
err = nf90_def_var(ncid,'salt',NCDF_FLOAT_PRECISION,f4_dims,salt_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,salt_id,long_name='salinity',units='PSU', &
FillValue=fv,missing_value=mv,valid_range=vr)
end if
if (save_s) then
fv = salt_missing
mv = salt_missing
vr(1) = 0.
vr(2) = 40.
err = nf90_def_var(ncid,'salt',NCDF_FLOAT_PRECISION,f4_dims,salt_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,salt_id,long_name='salinity',units='PSU', &
FillValue=fv,missing_value=mv,valid_range=vr)
end if
if (calc_temp .and. save_t) then
fv = temp_missing
mv = temp_missing
vr(1) = -2.
vr(2) = 40.
err = nf90_def_var(ncid,'temp',NCDF_FLOAT_PRECISION,f4_dims,temp_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,temp_id,long_name='temperature',units='degC',&
FillValue=fv,missing_value=mv,valid_range=vr)
end if
if (save_t) then
fv = temp_missing
mv = temp_missing
vr(1) = -2.
vr(2) = 40.
err = nf90_def_var(ncid,'temp',NCDF_FLOAT_PRECISION,f4_dims,temp_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,temp_id,long_name='temperature',units='degC',&
FillValue=fv,missing_value=mv,valid_range=vr)
end if
if (save_rho) then
fv = rho_missing
mv = rho_missing
vr(1) = 0.
vr(2) = 30.
err = nf90_def_var(ncid,'sigma_t',NCDF_FLOAT_PRECISION,f4_dims,sigma_t_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,sigma_t_id,long_name='sigma_t',units='kg/m3',&
FillValue=fv,missing_value=mv,valid_range=vr)
end if
if (save_rho) then
fv = rho_missing
mv = rho_missing
vr(1) = 0.
vr(2) = 30.
err = nf90_def_var(ncid,'sigma_t',NCDF_FLOAT_PRECISION,f4_dims,sigma_t_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,sigma_t_id,long_name='sigma_t',units='kg/m3',&
FillValue=fv,missing_value=mv,valid_range=vr)
end if
if (save_strho) then
if (save_rad) then
fv = rad_missing
mv = rad_missing
......@@ -328,7 +316,6 @@