!>@brief The module 'get_stochy_pattern_mod' contains the subroutines to retrieve the random pattern in the cubed-sphere grid
module get_stochy_pattern_mod
use kinddef
use spectral_transforms, only : len_trie_ls, &
len_trio_ls, ls_dim, stochy_la2ga, &
coslat_a, latg, levs, lonf, skeblevs,&
four_to_grid, spec_to_four, dezouv_stochy,dozeuv_stochy
use stochy_namelist_def, only : n_var_lndp, ntrunc, stochini,n_var_spp
use stochy_data_mod, only : gg_lats, gg_lons, inttyp, nskeb, nshum, nsppt, &
nocnsppt,nepbl,nlndp, &
rnlat, rpattern_sfc, rpattern_skeb, &
rpattern_shum, rpattern_sppt, rpattern_ocnsppt,&
rpattern_epbl1, rpattern_epbl2, skebu_save, &
nspp,rpattern_spp, &
skebv_save, skeb_vwts, skeb_vpts, wlon
use stochy_patterngenerator_mod, only: random_pattern, ndimspec, &
patterngenerator_advance
use stochy_internal_state_mod, only: stochy_internal_state
use mpi_wrapper, only : mp_reduce_sum,is_rootpe
use mersenne_twister, only: random_seed
implicit none
private
public get_random_pattern_vector,get_random_pattern_spp
public get_random_pattern_sfc,get_random_pattern_scalar
public write_stoch_restart_atm,write_stoch_restart_ocn
logical :: first_call=.true.
contains
!>@brief The subroutine 'get_random_pattern_sfc' converts spherical harmonics to the gaussian grid then interpolates to the target grid
!>@details This subroutine is for a 2-D (lat-lon) scalar field
subroutine get_random_pattern_sfc(rpattern,npatterns,&
gis_stochy,pattern_3d)
!\callgraph
! generate a random pattern for stochastic physics
implicit none
type(random_pattern), intent(inout) :: rpattern(npatterns)
type(stochy_internal_state), intent(in) :: gis_stochy
integer,intent(in):: npatterns
integer i,j,lat,n,k
real(kind=kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2d
! logical lprint
real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: workg
real (kind=kind_dbl_prec) glolal(lonf,gis_stochy%lats_node_a)
integer kmsk0(lonf,gis_stochy%lats_node_a)
real(kind=kind_dbl_prec),intent(out) :: pattern_3d(gis_stochy%nx,gis_stochy%ny,n_var_lndp)
real(kind=kind_dbl_prec) :: pattern_1d(gis_stochy%nx)
do k=1,n_var_lndp
kmsk0 = 0
glolal = 0.
do n=1,npatterns
call patterngenerator_advance(rpattern(n),k,.false.)
! if (is_rootpe()) print *, 'Random pattern for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(rpattern_sfc(n)%spec_o(:,:,k)), maxval(rpattern_sfc(n)%spec_o(:,:,k))
call scalarspect_to_gaugrid(rpattern(n),gis_stochy,wrk2d,k)
glolal = glolal + wrk2d(:,:,1)
enddo
allocate(workg(lonf,latg))
workg = 0.
do j=1,gis_stochy%lats_node_a
lat=gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+j)
do i=1,lonf
workg(i,lat) = glolal(i,j)
enddo
enddo
call mp_reduce_sum(workg,lonf,latg)
! if (is_rootpe()) print *, 'workg after mp_reduce_sum for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(workg), maxval(workg)
! interpolate to cube grid
do j=1,gis_stochy%ny
pattern_1d = 0
associate( tlats=>gis_stochy%parent_lats(1:gis_stochy%len(j),j),&
tlons=>gis_stochy%parent_lons(1:gis_stochy%len(j),j))
call stochy_la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,&
pattern_1d(1:gis_stochy%len(j)),gis_stochy%len(j),tlats,tlons)
pattern_3d(:,j,k)=pattern_1d(:)
end associate
enddo
! if (is_rootpe()) print *, '3D pattern for LNDP PERTS in get_random_pattern_fv3_sfc: k, min, max ',k,minval(pattern_3d(:,:,k)), maxval(pattern_3d(:,:,k))
deallocate(workg)
enddo ! loop over k, n_var_lndp
end subroutine get_random_pattern_sfc
!>@brief The subroutine 'get_random_pattern_fv3_vect' converts spherical harmonics to a vector on gaussian grid then interpolates to the target grid
!>@details This subroutine is for a 2-D (lat-lon) vector field
subroutine get_random_pattern_vector(rpattern,npatterns,&
gis_stochy,upattern_3d,vpattern_3d)
!\callgraph
! generate a random pattern for stochastic physics
implicit none
type(stochy_internal_state), intent(in) :: gis_stochy
type(random_pattern), intent(inout) :: rpattern(npatterns)
real(kind=kind_dbl_prec), dimension(len_trie_ls,2) :: vrtspec_e,divspec_e
real(kind=kind_dbl_prec), dimension(len_trio_ls,2) :: vrtspec_o,divspec_o
integer:: npatterns
real(kind=kind_dbl_prec) :: upattern_3d(gis_stochy%nx,gis_stochy%ny,levs)
real(kind=kind_dbl_prec) :: vpattern_3d(gis_stochy%nx,gis_stochy%ny,levs)
real(kind=kind_dbl_prec) :: pattern_1d(gis_stochy%nx)
integer i,j,lat,n,nn,k
real(kind_phys), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2du,wrk2dv
! logical lprint
real(kind_dbl_prec), allocatable, dimension(:,:) :: workgu,workgv
integer kmsk0(lonf,gis_stochy%lats_node_a)
kmsk0 = 0
allocate(workgu(lonf,latg))
allocate(workgv(lonf,latg))
divspec_e = 0; divspec_o = 0.
if (first_call) then
allocate(skebu_save(gis_stochy%nx,gis_stochy%ny,skeblevs))
allocate(skebv_save(gis_stochy%nx,gis_stochy%ny,skeblevs))
do k=2,skeblevs
workgu = 0.
workgv = 0.
do n=1,npatterns
if (.not. stochini) call patterngenerator_advance(rpattern(n),k,first_call)
! ke norm (convert streamfunction forcing to vorticity forcing)
do nn=1,2
vrtspec_e(:,nn) = gis_stochy%kenorm_e*rpattern(n)%spec_e(:,nn,k)
vrtspec_o(:,nn) = gis_stochy%kenorm_o*rpattern(n)%spec_o(:,nn,k)
enddo
! convert to winds
call vrtdivspect_to_uvgrid( divspec_e,divspec_o,vrtspec_e,vrtspec_o,&
wrk2du,wrk2dv, gis_stochy)
do i=1,lonf
do j=1,gis_stochy%lats_node_a
lat=gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+j)
workgu(i,lat) = workgu(i,lat) + wrk2du(i,j,1)
workgv(i,lat) = workgv(i,lat) + wrk2dv(i,j,1)
enddo
enddo
enddo
call mp_reduce_sum(workgu,lonf,latg)
call mp_reduce_sum(workgv,lonf,latg)
! interpolate to cube grid
do j=1,gis_stochy%ny
pattern_1d = 0
associate( tlats=>gis_stochy%parent_lats(1:gis_stochy%len(j),j),&
tlons=>gis_stochy%parent_lons(1:gis_stochy%len(j),j))
call stochy_la2ga(workgu,lonf,latg,gg_lons,gg_lats,wlon,rnlat,&
pattern_1d(1:gis_stochy%len(j)),gis_stochy%len(j),tlats,tlons)
skebu_save(:,j,k)=pattern_1d(:)
call stochy_la2ga(workgv,lonf,latg,gg_lons,gg_lats,wlon,rnlat,&
pattern_1d(1:gis_stochy%len(j)),gis_stochy%len(j),tlats,tlons)
skebv_save(:,j,k)=-1*pattern_1d(:)
end associate
enddo
enddo
endif
do k=1,skeblevs-1
skebu_save(:,:,k)=skebu_save(:,:,k+1)
skebv_save(:,:,k)=skebv_save(:,:,k+1)
do n=1,npatterns
rpattern(n)%spec_e(:,:,k)=rpattern(n)%spec_e(:,:,k+1)
rpattern(n)%spec_o(:,:,k)=rpattern(n)%spec_o(:,:,k+1)
enddo
enddo
! get pattern for last level
workgu = 0.
workgv = 0.
do n=1,npatterns
call patterngenerator_advance(rpattern(n),skeblevs,first_call)
! ke norm (convert streamfunction forcing to vorticity forcing)
divspec_e = 0; divspec_o = 0.
do nn=1,2
vrtspec_e(:,nn) = gis_stochy%kenorm_e*rpattern(n)%spec_e(:,nn,skeblevs)
vrtspec_o(:,nn) = gis_stochy%kenorm_o*rpattern(n)%spec_o(:,nn,skeblevs)
enddo
! convert to winds
call vrtdivspect_to_uvgrid(&
divspec_e,divspec_o,vrtspec_e,vrtspec_o,&
wrk2du,wrk2dv, gis_stochy)
do i=1,lonf
do j=1,gis_stochy%lats_node_a
lat=gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+j)
workgu(i,lat) = workgu(i,lat) + wrk2du(i,j,1)
workgv(i,lat) = workgv(i,lat) + wrk2dv(i,j,1)
enddo
enddo
enddo
call mp_reduce_sum(workgu,lonf,latg)
call mp_reduce_sum(workgv,lonf,latg)
! interpolate to cube grid
do j=1,gis_stochy%ny
pattern_1d = 0
associate( tlats=>gis_stochy%parent_lats(:,j),&
tlons=>gis_stochy%parent_lons(:,j))
call stochy_la2ga(workgu,lonf,latg,gg_lons,gg_lats,wlon,rnlat,&
pattern_1d(1:gis_stochy%len(j)),gis_stochy%len(j),tlats,tlons)
skebu_save(:,j,skeblevs)=pattern_1d(:)
call stochy_la2ga(workgv,lonf,latg,gg_lons,gg_lats,wlon,rnlat,&
pattern_1d(1:gis_stochy%len(j)),gis_stochy%len(j),tlats,tlons)
skebv_save(:,j,skeblevs)=-1*pattern_1d(:)
end associate
enddo
deallocate(workgu)
deallocate(workgv)
! interpolate in the vertical ! consider moving to cubed sphere side, more memory, but less interpolations
do k=1,levs
do j=1,gis_stochy%ny
upattern_3d(:,j,k) = skeb_vwts(k,1)*skebu_save(:,j,skeb_vpts(k,1))+skeb_vwts(k,2)*skebu_save(:,j,skeb_vpts(k,2))
vpattern_3d(:,j,k) = skeb_vwts(k,1)*skebv_save(:,j,skeb_vpts(k,1))+skeb_vwts(k,2)*skebv_save(:,j,skeb_vpts(k,2))
enddo
enddo
first_call=.false.
end subroutine get_random_pattern_vector
!>@brief The subroutine 'get_random_pattern_scalar' converts spherical harmonics to the gaussian grid then interpolates to the target grid
!>@details This subroutine is for a 2-D (lat-lon) scalar field
subroutine get_random_pattern_scalar(rpattern,npatterns,&
gis_stochy,pattern_2d)
! generate a random pattern for stochastic physics
implicit none
type(random_pattern), intent(inout) :: rpattern(npatterns)
type(stochy_internal_state) :: gis_stochy
integer,intent(in):: npatterns
integer i,j,lat,n
real(kind=kind_dbl_prec), dimension(lonf,gis_stochy%lats_node_a,1):: wrk2d
! logical lprint
real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: workg
real (kind=kind_dbl_prec) glolal(lonf,gis_stochy%lats_node_a)
integer kmsk0(lonf,gis_stochy%lats_node_a)
real(kind=kind_dbl_prec) :: pattern_2d(gis_stochy%nx,gis_stochy%ny)
real(kind=kind_dbl_prec) :: pattern_1d(gis_stochy%nx)
kmsk0 = 0
glolal = 0.
do n=1,npatterns
call patterngenerator_advance(rpattern(n),1,.false.)
call scalarspect_to_gaugrid(rpattern(n),gis_stochy, &
wrk2d,1)
glolal = glolal + wrk2d(:,:,1)
enddo
allocate(workg(lonf,latg))
workg = 0.
do j=1,gis_stochy%lats_node_a
lat=gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+j)
do i=1,lonf
workg(i,lat) = glolal(i,j)
enddo
enddo
call mp_reduce_sum(workg,lonf,latg)
! interpolate to cube grid
do j=1,gis_stochy%ny
pattern_1d = 0
associate( tlats=>gis_stochy%parent_lats(1:gis_stochy%len(j),j),&
tlons=>gis_stochy%parent_lons(1:gis_stochy%len(j),j))
call stochy_la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,&
pattern_1d(1:gis_stochy%len(j)),gis_stochy%len(j),tlats,tlons)
pattern_2d(:,j)=pattern_1d(:)
end associate
enddo
deallocate(workg)
end subroutine get_random_pattern_scalar
!>@brief The subroutine 'get_random_pattern_spp' converts spherical harmonics
!to the gaussian grid then interpolates to the target grid
!>@details This subroutine is for a 2-D (lat-lon) scalar field
subroutine get_random_pattern_spp(rpattern,npatterns,&
gis_stochy,pattern_3d)
! generate a random pattern for stochastic physics
implicit none
type(random_pattern), intent(inout) :: rpattern(npatterns)
type(stochy_internal_state) :: gis_stochy
integer,intent(in):: npatterns
integer i,j,lat,n
! logical lprint
real(kind=kind_dbl_prec), allocatable, dimension(:,:) :: workg
real (kind=kind_dbl_prec) glolal(lonf,gis_stochy%lats_node_a)
integer kmsk0(lonf,gis_stochy%lats_node_a)
real(kind=kind_dbl_prec) :: pattern_3d(gis_stochy%nx,gis_stochy%ny,npatterns)
real(kind=kind_dbl_prec) :: pattern_1d(gis_stochy%nx)
allocate(workg(lonf,latg))
do n=1,npatterns
kmsk0 = 0
glolal = 0.
call patterngenerator_advance(rpattern(n),1,.false.)
call scalarspect_to_gaugrid(rpattern(n),gis_stochy, &
glolal,1)
workg = 0.
do j=1,gis_stochy%lats_node_a
lat=gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+j)
do i=1,lonf
workg(i,lat) = glolal(i,j)
enddo
enddo
call mp_reduce_sum(workg,lonf,latg)
! interpolate to cube grid
do j=1,gis_stochy%ny
pattern_1d = 0
associate( tlats=>gis_stochy%parent_lats(1:gis_stochy%len(j),j),&
tlons=>gis_stochy%parent_lons(1:gis_stochy%len(j),j))
call stochy_la2ga(workg,lonf,latg,gg_lons,gg_lats,wlon,rnlat,&
pattern_1d(1:gis_stochy%len(j)),gis_stochy%len(j),tlats,tlons)
pattern_3d(:,j,n)=pattern_1d(:)
end associate
enddo
enddo
deallocate(workg)
end subroutine get_random_pattern_spp
!>@brief The subroutine 'scalarspect_to_gaugrid' converts scalar spherical harmonics to a scalar on a gaussian grid
!>@details This subroutine is for a 2-D (lat-lon) scalar field
subroutine scalarspect_to_gaugrid(rpattern,gis_stochy,datag,n)
!\callgraph
implicit none
type(random_pattern), intent(in) :: rpattern
type(stochy_internal_state), intent(in) :: gis_stochy
integer , intent(in) :: n
real(kind=kind_dbl_prec), intent(out) :: datag(lonf,gis_stochy%lats_node_a)
! local vars
real(kind=kind_dbl_prec) for_gr_a_1(gis_stochy%lon_dim_a,1,gis_stochy%lats_dim_a)
real(kind=kind_dbl_prec) for_gr_a_2(lonf,1,gis_stochy%lats_dim_a)
integer i,k
integer lan,lat
call spec_to_four(rpattern%spec_e(:,:,n), rpattern%spec_o(:,:,n), &
gis_stochy%plnev_a,gis_stochy%plnod_a,&
gis_stochy%ls_node, &
gis_stochy%lats_dim_a,for_gr_a_1,&
gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,&
gis_stochy%lats_nodes_a,gis_stochy%global_lats_a,&
gis_stochy%lats_node_a,gis_stochy%ipt_lats_node_a,1)
do lan=1,gis_stochy%lats_node_a
lat = gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+lan)
call four_to_grid(for_gr_a_1(:,:,lan),for_gr_a_2(:,:,lan),&
gis_stochy%lon_dim_a,1)
enddo
datag = 0.
do lan=1,gis_stochy%lats_node_a
lat = gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+lan)
do i=1,lonf
datag(i,lan) = for_gr_a_2(i,1,lan)
enddo
enddo
return
end subroutine scalarspect_to_gaugrid
!>@brief The subroutine 'write_patterns' writes out a single pattern and the seed associated with the random number sequence in netcdf
!>@brief The subroutine 'write_stoch_restart_atm' writes out the speherical harmonics to a file, controlled by restart_interval
!>@details Only the active patterns are written out
subroutine write_stoch_restart_atm(sfile)
!\callgraph
use netcdf
use stochy_namelist_def, only : do_sppt,do_shum,do_skeb,lndp_type,do_spp
implicit none
character(len=*) :: sfile
integer :: stochlun,k,n,isize,ierr
integer :: ncid,varid1a,varid1b,varid2a,varid2b,varid3a,varid3b,varid4a,varid4b,varid5a,varid5b
integer :: seed_dim_id,spec_dim_id,zt_dim_id,ztsfc_dim_id,np_dim_id,npsfc_dim_id
integer :: ztspp_dim_id,npspp_dim_id
include 'netcdf.inc'
if ( ( .NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (lndp_type==0 ) .AND. (.NOT. do_spp)) return
stochlun=99
if (is_rootpe()) then
if (nsppt > 0 .OR. nshum > 0 .OR. nskeb > 0 .OR. nlndp>0 .OR. nspp>0 ) then
ierr=nf90_create(trim(sfile),cmode=NF90_CLOBBER,ncid=ncid)
ierr=NF90_PUT_ATT(ncid,NF_GLOBAL,"ntrunc",ntrunc)
call random_seed(size=isize) ! get seed size
ierr=NF90_DEF_DIM(ncid,"len_seed",isize,seed_dim_id)
ierr=NF90_PUT_ATT(ncid,seed_dim_id,"long_name","length of random seed")
ierr=NF90_DEF_DIM(ncid,"num_patterns",NF_UNLIMITED,np_dim_id) ! should be 5
ierr=NF90_PUT_ATT(ncid,np_dim_id,"long_name","number of random patterns (max of 5)")
if (lndp_type .NE. 0) then
ierr=NF90_DEF_DIM(ncid,"num_patterns_sfc",nlndp,npsfc_dim_id) ! should be 5
ierr=NF90_PUT_ATT(ncid,npsfc_dim_id,"long_name","number of random patterns for surface)")
ierr=NF90_DEF_DIM(ncid,"n_var_lndp",n_var_lndp,ztsfc_dim_id)
ierr=NF90_PUT_ATT(ncid,ztsfc_dim_id,"long_name","number of sfc perturbation types")
endif
if (nspp .GT. 0) then
ierr=NF90_DEF_DIM(ncid,"num_patterns_spp",nspp,npspp_dim_id) ! should be 5
ierr=NF90_PUT_ATT(ncid,npspp_dim_id,"long_name","number of random patterns for spp)")
ierr=NF90_DEF_DIM(ncid,"n_var_spp",n_var_spp,ztspp_dim_id)
ierr=NF90_PUT_ATT(ncid,ztspp_dim_id,"long_name","number of spp perturbation types")
endif
ierr=NF90_DEF_DIM(ncid,"ndimspecx2",2*ndimspec,spec_dim_id)
ierr=NF90_PUT_ATT(ncid,spec_dim_id,"long_name","number of spectral cofficients")
if (do_sppt) then
ierr=NF90_DEF_VAR(ncid,"sppt_seed",NF90_DOUBLE,(/seed_dim_id, np_dim_id/), varid1a)
ierr=NF90_PUT_ATT(ncid,varid1a,"long_name","random number seed for SPPT")
ierr=NF90_DEF_VAR(ncid,"sppt_spec",NF90_DOUBLE,(/spec_dim_id, np_dim_id/), varid1b)
ierr=NF90_PUT_ATT(ncid,varid1b,"long_name","spectral cofficients SPPT")
endif
if (do_shum) then
ierr=NF90_DEF_VAR(ncid,"shum_seed",NF90_DOUBLE,(/seed_dim_id, np_dim_id/), varid2a)
ierr=NF90_PUT_ATT(ncid,varid2a,"long_name","random number seed for SHUM")
ierr=NF90_DEF_VAR(ncid,"shum_spec",NF90_DOUBLE,(/spec_dim_id, np_dim_id/), varid2b)
ierr=NF90_PUT_ATT(ncid,varid2b,"long_name","spectral cofficients SHUM")
endif
if (do_skeb) then
ierr=NF90_DEF_DIM(ncid,"skeblevs",skeblevs,zt_dim_id)
ierr=NF90_PUT_ATT(ncid,zt_dim_id,"long_name","number of vertical levels for SKEB")
ierr=NF90_DEF_VAR(ncid,"skeb_seed",NF90_DOUBLE,(/seed_dim_id, zt_dim_id,np_dim_id/), varid3a)
ierr=NF90_PUT_ATT(ncid,varid3a,"long_name","random number seed for SKEB")
ierr=NF90_DEF_VAR(ncid,"skeb_spec",NF90_DOUBLE,(/spec_dim_id, zt_dim_id,np_dim_id/), varid3b)
ierr=NF90_PUT_ATT(ncid,varid3b,"long_name","spectral cofficients SKEB")
endif
if (nlndp>0) then
ierr=NF90_DEF_VAR(ncid,"sfcpert_seed",NF90_DOUBLE,(/seed_dim_id, ztsfc_dim_id, npsfc_dim_id/), varid4a)
ierr=NF90_PUT_ATT(ncid,varid4a,"long_name","random number seed for SHUM")
ierr=NF90_DEF_VAR(ncid,"sfcpert_spec",NF90_DOUBLE,(/spec_dim_id, ztsfc_dim_id, npsfc_dim_id/), varid4b)
ierr=NF90_PUT_ATT(ncid,varid4b,"long_name","spectral cofficients SHUM")
endif
if (nspp>0) then
ierr=NF90_DEF_VAR(ncid,"spp_seed",NF90_DOUBLE,(/seed_dim_id, ztspp_dim_id, npspp_dim_id/), varid5a)
ierr=NF90_PUT_ATT(ncid,varid5a,"long_name","random number seed for SPP")
ierr=NF90_DEF_VAR(ncid,"spp_spec",NF90_DOUBLE,(/spec_dim_id, ztspp_dim_id, npspp_dim_id/), varid5b)
ierr=NF90_PUT_ATT(ncid,varid5b,"long_name","spectral cofficients SPP")
endif
ierr=NF90_ENDDEF(ncid)
if (ierr .NE. 0) then
write(0,*) 'error creating stochastic restart file'
return
end if
endif
endif
if (nsppt > 0) then
do n=1,nsppt
call write_pattern(rpattern_sppt(n),ncid,1,n,varid1a,varid1b,.false.,ierr)
enddo
endif
if (nshum > 0) then
do n=1,nshum
call write_pattern(rpattern_shum(n),ncid,1,n,varid2a,varid2b,.false.,ierr)
enddo
endif
if (nskeb > 0) then
do n=1,nskeb
do k=1,skeblevs
call write_pattern(rpattern_skeb(n),ncid,k,n,varid3a,varid3b,.true.,ierr)
enddo
enddo
endif
if (lndp_type .NE. 0 .AND. nlndp>0) then
do n=1,nlndp
do k=1,n_var_lndp
call write_pattern(rpattern_sfc(n),ncid,k,n,varid4a,varid4b,.true.,ierr)
enddo
enddo
endif
if (nspp > 0) then
do n=1,nspp
call write_pattern(rpattern_spp(n),ncid,1,n,varid5a,varid5b,.true.,ierr)
enddo
endif
if (is_rootpe() ) then
ierr=NF90_CLOSE(ncid)
if (ierr .NE. 0) then
write(0,*) 'error writing patterns and closing file'
return
endif
endif
end subroutine write_stoch_restart_atm
!>@brief The subroutine 'write_stoch_restart_ocn' writes out the speherical harmonics to a file, controlled by restart_interval
!>@details Only the active patterns are written out
subroutine write_stoch_restart_ocn(sfile)
!\callgraph
use netcdf
use stochy_namelist_def, only : do_ocnsppt,pert_epbl
implicit none
character(len=*) :: sfile
integer :: stochlun,k,n,isize,ierr
integer :: ncid,varid1a,varid1b,varid2a,varid2b,varid3a,varid3b
integer :: seed_dim_id,spec_dim_id,np_dim_id
include 'netcdf.inc'
if ( ( .NOT. do_ocnsppt) .AND. (.NOT. pert_epbl) ) return
stochlun=99
if (is_rootpe()) then
ierr=nf90_create(trim(sfile),cmode=NF90_CLOBBER,ncid=ncid)
ierr=NF90_PUT_ATT(ncid,NF_GLOBAL,"ntrunc",ntrunc)
call random_seed(size=isize) ! get seed size
ierr=NF90_DEF_DIM(ncid,"len_seed",isize,seed_dim_id)
ierr=NF90_PUT_ATT(ncid,seed_dim_id,"long_name","length of random seed")
ierr=NF90_DEF_DIM(ncid,"num_patterns",NF_UNLIMITED,np_dim_id) ! should be 5
ierr=NF90_PUT_ATT(ncid,np_dim_id,"long_name","number of random patterns (max of 5)")
ierr=NF90_DEF_DIM(ncid,"ndimspecx2",2*ndimspec,spec_dim_id)
ierr=NF90_PUT_ATT(ncid,spec_dim_id,"long_name","number of spectral cofficients")
if (do_ocnsppt) then
ierr=NF90_DEF_VAR(ncid,"ocnsppt_seed",NF90_DOUBLE,(/seed_dim_id, np_dim_id/), varid1a)
ierr=NF90_PUT_ATT(ncid,varid1a,"long_name","random number seed for SPPT")
ierr=NF90_DEF_VAR(ncid,"ocnsppt_spec",NF90_DOUBLE,(/spec_dim_id, np_dim_id/), varid1b)
ierr=NF90_PUT_ATT(ncid,varid1b,"long_name","spectral cofficients SPPT")
endif
if (pert_epbl) then
ierr=NF90_DEF_VAR(ncid,"epbl1_seed",NF90_DOUBLE,(/seed_dim_id, np_dim_id/), varid2a)
ierr=NF90_PUT_ATT(ncid,varid2a,"long_name","random number seed for EPBL1")
ierr=NF90_DEF_VAR(ncid,"epbl1_spec",NF90_DOUBLE,(/spec_dim_id, np_dim_id/), varid2b)
ierr=NF90_PUT_ATT(ncid,varid2b,"long_name","spectral cofficients EPBL1")
ierr=NF90_DEF_VAR(ncid,"epbl2_seed",NF90_DOUBLE,(/seed_dim_id, np_dim_id/), varid3a)
ierr=NF90_PUT_ATT(ncid,varid3a,"long_name","random number seed for EPBL2")
ierr=NF90_DEF_VAR(ncid,"epbl2_spec",NF90_DOUBLE,(/spec_dim_id, np_dim_id/), varid3b)
ierr=NF90_PUT_ATT(ncid,varid3b,"long_name","spectral cofficients EPBL2")
endif
ierr=NF90_ENDDEF(ncid)
if (ierr .NE. 0) then
write(0,*) 'error creating stochastic restart file'
return
end if
endif
if (nocnsppt > 0) then
do n=1,nocnsppt
call write_pattern(rpattern_ocnsppt(n),ncid,1,n,varid1a,varid1b,.false.,ierr)
enddo
endif
if (nepbl > 0) then
do n=1,nepbl
call write_pattern(rpattern_epbl1(n),ncid,1,n,varid2a,varid2b,.false.,ierr)
call write_pattern(rpattern_epbl2(n),ncid,1,n,varid3a,varid3b,.false.,ierr)
enddo
endif
if (is_rootpe() ) then
ierr=NF90_CLOSE(ncid)
if (ierr .NE. 0) then
write(0,*) 'error writing patterns and closing file'
return
endif
endif
end subroutine write_stoch_restart_ocn
!>@brief The subroutine 'write_patterns' writes out a single pattern and the seed associated with the random number sequence
!>@details Spherical harminoncs are stored with trianglular truncation
subroutine write_pattern(rpattern,outlun,lev,np,varid1,varid2,slice_of_3d,iret)
!\callgraph
use netcdf
implicit none
type(random_pattern), intent(inout) :: rpattern
integer, intent(in) :: outlun,lev
integer, intent(in) :: np,varid1,varid2
logical, intent(in) :: slice_of_3d
integer, intent(out) :: iret
real(kind_phys), allocatable :: pattern2d(:)
integer nm,nn,arrlen,isize,ierr
integer,allocatable :: isave(:)
include 'netcdf.inc'
arrlen=2*ndimspec
iret=0
allocate(pattern2d(arrlen))
pattern2d=0.0
! fill in apprpriate pieces of array
do nn=1,len_trie_ls
nm = rpattern%idx_e(nn)
if (nm == 0) cycle
pattern2d(nm) = rpattern%spec_e(nn,1,lev)
pattern2d(ndimspec+nm) = rpattern%spec_e(nn,2,lev)
enddo
do nn=1,len_trio_ls
nm = rpattern%idx_o(nn)
if (nm == 0) cycle
pattern2d(nm) = rpattern%spec_o(nn,1,lev)
pattern2d(ndimspec+nm) = rpattern%spec_o(nn,2,lev)
enddo
call mp_reduce_sum(pattern2d,arrlen)
! write only on root process
if (is_rootpe()) then
print*,'writing out random pattern (min/max/size)',&
minval(pattern2d),maxval(pattern2d),size(pattern2d)
call random_seed(size=isize) ! get seed size
allocate(isave(isize)) ! get seed
call random_seed(get=isave,stat=rpattern%rstate) ! write seed
ierr=NF90_PUT_VAR(outlun,varid1,isave,(/1,np/))
if (slice_of_3d) then
ierr=NF90_PUT_VAR(outlun,varid2,pattern2d,(/1,lev,np/))
else
ierr=NF90_PUT_VAR(outlun,varid2,pattern2d,(/1,np/))
endif
if (ierr .NE. 0) then
write(0,*) 'error writing to stochastic restart file'
iret = ierr
return
end if
endif
deallocate(pattern2d)
end subroutine write_pattern
!>@brief The subroutine 'vrtdivspect_to_uvgrid' converts vorticty and divergence spherical harmonics to
! zonal and meridional winds on the gaussian grid
!>@details This subroutine is for a 2-D (lat-lon) vector field
subroutine vrtdivspect_to_uvgrid(&
trie_di,trio_di,trie_ze,trio_ze,&
uug,vvg, gis_stochy)
!\callgraph
implicit none
type(stochy_internal_state), intent(in) :: gis_stochy
real(kind=kind_dbl_prec), intent(in) :: trie_di(len_trie_ls,2)
real(kind=kind_dbl_prec), intent(in) :: trio_di(len_trio_ls,2)
real(kind=kind_dbl_prec), intent(in) :: trie_ze(len_trie_ls,2)
real(kind=kind_dbl_prec), intent(in) :: trio_ze(len_trio_ls,2)
real(kind=kind_phys), intent(out) :: uug(lonf,gis_stochy%lats_node_a)
real(kind=kind_phys), intent(out) :: vvg(lonf,gis_stochy%lats_node_a)
! local vars
real(kind=kind_dbl_prec) trie_ls(len_trie_ls,2,2)
real(kind=kind_dbl_prec) trio_ls(len_trio_ls,2,2)
real(kind=kind_dbl_prec) for_gr_a_1(gis_stochy%lon_dim_a,2,gis_stochy%lats_dim_a)
real(kind=kind_dbl_prec) for_gr_a_2(lonf,2,gis_stochy%lats_dim_a)
integer i,k
integer lan,lat
real (kind=kind_phys) tx1
call dezouv_stochy(trie_di(:,:), trio_ze(:,:), &
trie_ls(:,:,1), trio_ls(:,:,2), gis_stochy%epsedn,gis_stochy%epsodn, &
gis_stochy%snnp1ev,gis_stochy%snnp1od,gis_stochy%ls_node)
call dozeuv_stochy(trio_di(:,:), trie_ze(:,:), &
trio_ls(:,:,1), trie_ls(:,:,2), gis_stochy%epsedn,gis_stochy%epsodn, &
gis_stochy%snnp1ev,gis_stochy%snnp1od,gis_stochy%ls_node)
call spec_to_four(trie_ls, trio_ls, &
gis_stochy%plnev_a,gis_stochy%plnod_a,&
gis_stochy%ls_node,&
gis_stochy%lats_dim_a,for_gr_a_1,&
gis_stochy%ls_nodes,gis_stochy%max_ls_nodes,&
gis_stochy%lats_nodes_a,gis_stochy%global_lats_a,&
gis_stochy%lats_node_a,gis_stochy%ipt_lats_node_a,2)
do lan=1,gis_stochy%lats_node_a
lat = gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+lan)
call four_to_grid(for_gr_a_1(:,:,lan),for_gr_a_2(:,:,lan),&
gis_stochy%lon_dim_a,2)
enddo
uug = 0.; vvg = 0.
do lan=1,gis_stochy%lats_node_a
lat = gis_stochy%global_lats_a(gis_stochy%ipt_lats_node_a-1+lan)
tx1 = 1. / coslat_a(lat)
do i=1,lonf
uug(i,lan) = for_gr_a_2(i,1,lan) * tx1
vvg(i,lan) = for_gr_a_2(i,2,lan) * tx1
enddo
enddo
return
end subroutine vrtdivspect_to_uvgrid
end module get_stochy_pattern_mod