diff --git a/Makefile b/Makefile
index 5a4f8fd..1226065 100644
--- a/Makefile
+++ b/Makefile
@@ -49,7 +49,7 @@ KNDSTR = \
set_knd.o
OBJSTR = \
filename_mod.o\
- da_params.o\
+ da_params.o\
drv_str.o\
cns_str.o\
obs_str.o\
@@ -100,11 +100,13 @@ PHYSOBS = \
OBJS = \
def_nml.o\
+ def_nml_multi.o\
def_grd.o\
sav_itr.o\
rdeofs_chl.o\
rdeofs_n3n.o\
rdeofs_o2o.o\
+ rdeofs_multi.o\
rdrcorr.o\
mean_rdr.o\
netcdf_err.o\
@@ -129,6 +131,7 @@ OBJS = \
obs_arg_ad.o\
veof_chl_ad.o\
veof_nut_ad.o\
+ veof_multiv_ad.o\
ver_hor_chl_ad.o\
ver_hor_nut_ad.o\
rcfl_x_ad.o\
@@ -146,8 +149,11 @@ OBJS = \
readChlStat.o\
readNutStat.o\
readNutCov.o\
+ readChlNutCov.o\
wrt_chl_stat.o\
+ wrt_upd_nut.o\
wrt_nut_stat.o\
+ cp_nut_stat.o\
costf.o\
obs_sat.o\
bio_conv.o\
diff --git a/bio_conv.f90 b/bio_conv.f90
index cffb5fb..4f59adf 100644
--- a/bio_conv.f90
+++ b/bio_conv.f90
@@ -29,16 +29,22 @@ subroutine bio_conv
!-----------------------------------------------------------------------
use grd_str
+ use eof_str
use bio_str
+ use drv_str
implicit none
- INTEGER(i4) :: i, j, k, l
+ INTEGER(i4) :: i, j, k, l, my_km
grd%chl(:,:,:) = 0.0
+ my_km = grd%km
+ if(drv%multiv.eq.1) &
+ my_km = ros%kmchl
+
do l = 1,bio%nphy
- do k = 1,grd%km
+ do k = 1,my_km
do j = 1,grd%jm
do i = 1,grd%im
grd%chl(i,j,k) = grd%chl(i,j,k) + bio%phy(i,j,k,l,1)
diff --git a/bio_conv_ad.f90 b/bio_conv_ad.f90
index 4378c9e..6d2582e 100644
--- a/bio_conv_ad.f90
+++ b/bio_conv_ad.f90
@@ -29,16 +29,22 @@ subroutine bio_conv_ad
!-----------------------------------------------------------------------
use grd_str
+ use eof_str
use bio_str
+ use drv_str
implicit none
- INTEGER(i4) :: i, j, k, l
+ INTEGER(i4) :: i, j, k, l, my_km
bio%phy_ad(:,:,:,:,:) = 0.0
+ my_km = grd%km
+ if(drv%multiv .eq. 1) &
+ my_km = ros%kmchl
+
do l = 1,bio%nphy
- do k = 1,grd%km
+ do k = 1,my_km
do j = 1,grd%jm
do i = 1,grd%im
bio%phy_ad(i,j,k,l,1) = bio%phy_ad(i,j,k,l,1) + grd%chl_ad(i,j,k)
diff --git a/bio_mod.f90 b/bio_mod.f90
index 7d10468..95b3f48 100644
--- a/bio_mod.f90
+++ b/bio_mod.f90
@@ -31,16 +31,21 @@ subroutine bio_mod
use grd_str
use bio_str
use eof_str
+ use drv_str
IMPLICIT NONE
- INTEGER(i4) :: m, l, k,j ,i
+ INTEGER(i4) :: m, l, k,j ,i, my_km
bio%phy(:,:,:,:,:) = 0.0
+ my_km = grd%km
+ if(drv%multiv .eq. 1) &
+ my_km = ros%kmchl
+
do m=1,bio%ncmp
do l=1,bio%nphy
- do k=1,grd%km
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
bio%phy(i,j,k,l,m)=bio%cquot(i,j,k,l,m)*bio%pquot(i,j,k,l)*grd%chl(i,j,k)
diff --git a/bio_mod_ad.f90 b/bio_mod_ad.f90
index 1ab9378..d714ae0 100644
--- a/bio_mod_ad.f90
+++ b/bio_mod_ad.f90
@@ -31,16 +31,21 @@ subroutine bio_mod_ad
use grd_str
use bio_str
use eof_str
+ use drv_str
IMPLICIT NONE
- INTEGER(i4) :: m, l, k, j, i
+ INTEGER(i4) :: m, l, k, j, i, my_km
grd%chl_ad(:,:,:) = 0.0
+ my_km = grd%km
+ if(drv%multiv .eq. 1) &
+ my_km = ros%kmchl
+
do m=1,bio%ncmp
do l=1,bio%nphy
- do k=1,grd%km
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
grd%chl_ad(i,j,k) = grd%chl_ad(i,j,k) + bio%cquot(i,j,k,l,m) * bio%pquot(i,j,k,l) * bio%phy_ad(i,j,k,l,m)
diff --git a/bio_str.f90 b/bio_str.f90
index 9d1b642..724ac21 100644
--- a/bio_str.f90
+++ b/bio_str.f90
@@ -42,7 +42,9 @@ MODULE bio_str
REAL(r8), POINTER :: phy_ad(:,:,:,:,:) ! biogeochemical adjoint variables
REAL(r8), POINTER :: InitialChl(:,:,:) ! initial amount of chlorophyll
REAL(r8), POINTER :: InitialNut(:,:,:,:) ! initial amount of nutrients
- REAL(r8), POINTER :: covn3n_n1p(:,:,:) ! initial amount of nutrients
+ REAL(r8), POINTER :: covn3n_n1p(:,:,:) ! covariance n3n n1p
+ REAL(r8), POINTER :: covn3n_chl(:,:,:) ! covariance n3n chl
+ REAL(r8), POINTER :: covn1p_chl(:,:,:) ! covariance n3n chl
INTEGER :: nphy ! number of phytoplankton types
INTEGER :: ncmp ! No. of phytoplankton components
diff --git a/clean_mem.f90 b/clean_mem.f90
index bdbf1ae..a05255b 100644
--- a/clean_mem.f90
+++ b/clean_mem.f90
@@ -73,11 +73,16 @@ subroutine clean_mem
DEALLOCATE ( rcf%sc)
DEALLOCATE(SendCountX3D, SendDisplX3D)
+ DEALLOCATE(SendCountX3D_chl, SendDisplX3D_chl)
DEALLOCATE(RecCountX3D, RecDisplX3D)
+ DEALLOCATE(RecCountX3D_chl, RecDisplX3D_chl)
DEALLOCATE(ChlExtended)
+ DEALLOCATE(ChlExtended_3d,N3nExtended_3d,O2oExtended_3d)
DEALLOCATE(SendBottom, RecTop)
DEALLOCATE(SendTop, RecBottom)
+ DEALLOCATE(SendTop_2d, RecBottom_2d)
+ DEALLOCATE(SendBottom_2d, RecTop_2d)
if(MyId .eq. 0) then
write(*,*) ' ALL MEMORY CLEAN'
diff --git a/cnv_inn.f90 b/cnv_inn.f90
index 7f098f5..ca686f0 100644
--- a/cnv_inn.f90
+++ b/cnv_inn.f90
@@ -43,21 +43,28 @@ subroutine cnv_inn
! Convert the control vector to v
call cnv_ctv
- if(drv%chl_assim .eq. 1) then
- call ver_hor_chl
- endif
- if(drv%nut .eq. 1) then
- if(bio%N3n .eq. 1) then
- call ver_hor_nut(grd%n3n, grd%n3n_ad,'N')
+ if (drv%multiv .eq. 0) then
+ if(drv%chl_assim .eq. 1) then
+ call ver_hor_chl
endif
- if(bio%O2o .eq. 1) then
- call ver_hor_nut(grd%o2o, grd%o2o_ad,'O')
+ if(drv%nut .eq. 1) then
+ if(bio%N3n .eq. 1) then
+ call ver_hor_nut(grd%n3n, grd%n3n_ad,'N')
+ endif
+ if(bio%O2o .eq. 1) then
+ call ver_hor_nut(grd%o2o, grd%o2o_ad,'O')
+ endif
endif
endif
+ if (drv%multiv .eq. 1) then
+ call ver_hor_chl
+ call ver_hor_nut(grd%n3n, grd%n3n_ad,'N')
+ endif
+
! ---
! Apply biological repartition of the chlorophyll
- if(drv%chl_assim .eq. 1) &
+ if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) &
call bio_mod
end subroutine cnv_inn
diff --git a/costf.f90 b/costf.f90
index 332435c..41570c1 100644
--- a/costf.f90
+++ b/costf.f90
@@ -55,22 +55,28 @@ subroutine costf
call cnv_ctv
! --------
- ! Control to physical space
- if(drv%chl_assim .eq. 1) then
- call ver_hor_chl
- endif
- if(drv%nut .eq. 1) then
- if(bio%N3n .eq. 1) then
- call ver_hor_nut(grd%n3n, grd%n3n_ad, 'N')
+ ! Control to physical space
+ if (drv%multiv.eq.0) then
+ if(drv%chl_assim .eq. 1) then
+ call ver_hor_chl
endif
- if(bio%O2o .eq. 1) then
- call ver_hor_nut(grd%o2o, grd%o2o_ad, 'O')
+ if(drv%nut .eq. 1) then
+ if(bio%N3n .eq. 1) then
+ call ver_hor_nut(grd%n3n, grd%n3n_ad, 'N')
+ endif
+ if(bio%O2o .eq. 1) then
+ call ver_hor_nut(grd%o2o, grd%o2o_ad, 'O')
+ endif
endif
+
+ else if(drv%multiv .eq. 1) then
+ call ver_hor_chl
+ call ver_hor_nut(grd%n3n, grd%n3n_ad, 'N')
endif
! ---
! Apply biological repartition of the chlorophyll
- if(drv%chl_assim .eq. 1) &
+ if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) &
call bio_mod
! --------
@@ -109,23 +115,29 @@ subroutine costf
! ---
! Apply biological repartition of the chlorophyll
- if(drv%chl_assim .eq. 1) &
+ if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) &
call bio_mod_ad
! --------
- ! Control to physical space
- if(drv%chl_assim .eq. 1) then
- call ver_hor_chl_ad
- endif
- if(drv%nut .eq. 1) then
- if(bio%N3n .eq. 1) then
- call ver_hor_nut_ad(grd%n3n, grd%n3n_ad, 'N')
+ ! Control to physical space
+ if (drv%multiv .eq. 0) then
+ if(drv%chl_assim .eq. 1) then
+ call ver_hor_chl_ad
endif
- if(bio%O2o .eq. 1) then
- call ver_hor_nut_ad(grd%o2o, grd%o2o_ad, 'O')
+ if(drv%nut .eq. 1) then
+ if(bio%N3n .eq. 1) then
+ call ver_hor_nut_ad(grd%n3n, grd%n3n_ad, 'N')
+ endif
+ if(bio%O2o .eq. 1) then
+ call ver_hor_nut_ad(grd%o2o, grd%o2o_ad, 'O')
+ endif
endif
+
+ else if(drv%multiv .eq. 1) then
+ call ver_hor_chl_ad
+ call ver_hor_nut_ad(grd%n3n, grd%n3n_ad, 'N')
+ call veof_multiv_ad
endif
-
! write(*,*) 'COSTF sum(ro_ad) = ' , sum(grd%ro_ad)
! --------
! Convert the control vector
diff --git a/cp_nut_stat.f90 b/cp_nut_stat.f90
new file mode 100644
index 0000000..23e5c2b
--- /dev/null
+++ b/cp_nut_stat.f90
@@ -0,0 +1,127 @@
+subroutine cp_nut_stat
+
+ use set_knd
+ use grd_str
+ use drv_str
+ use mpi_str
+ use bio_str
+ use pnetcdf
+ use da_params
+
+ implicit none
+
+ INTEGER(i4) :: ncid, ierr, i, j, k, l
+ INTEGER(i4) :: idP, iVar
+ INTEGER(I4) :: xid,yid,depid,timeId, idTim
+ INTEGER :: system, SysErr
+
+ INTEGER(kind=MPI_OFFSET_KIND) :: global_im, global_jm, global_km, MyTime
+ INTEGER(KIND=MPI_OFFSET_KIND) :: MyCountSingle(1), MyStartSingle(1)
+ CHARACTER(LEN=46) :: BioRestart
+ CHARACTER(LEN=47) :: BioRestartLong
+ CHARACTER(LEN=6) :: MyVarName
+ ! LOGICAL, ALLOCATABLE :: MyConditions(:,:,:,:)
+
+ ! bug fix Intel 2018
+ real(r4), allocatable, dimension(:,:,:,:) :: DumpBio
+ integer(KIND=MPI_OFFSET_KIND) :: MyStart_4d(4), MyCount_4d(4)
+
+ real(r8) :: TimeArr(1)
+
+
+ MyStart_4d(1:3) = MyStart(:)
+ MyStart_4d(4) = 1
+ MyCount_4d(1:3) = MyCount(:)
+ MyCount_4d(4) = 1
+
+
+ ALLOCATE(DumpBio(grd%im,grd%jm,grd%km,1)); DumpBio(:,:,:,:) = 1.e20
+ ! ALLOCATE(MyConditions(grd%im,grd%jm,grd%km,bio%nphy))
+
+ if(MyId .eq. 0) then
+ write(drv%dia,*) 'writing nut structure (only copy from RSTbefore)'
+ write(*,*) 'writing nut structure (only copy from RSTbefore)'
+ endif
+
+ global_im = GlobalRow
+ global_jm = GlobalCol
+ global_km = grd%km
+ MyTime = 1
+
+ MyCountSingle(1) = 1
+ MyStartSingle(1) = 1
+ TimeArr(1) = DA_JulianDate
+
+
+ do l=1,NNutVar
+ iVar = NPhytoVar + l
+
+ if(iVar .gt. NBioVar) then
+ if(MyId .eq. 0) &
+ write(*,*) "Warning: Reading a variable not in the DA_VarList!"
+ endif
+
+ BioRestart = 'DA__FREQ_1/RST_after.'//ShortDate//'.'//DA_VarList(iVar)//'.nc'
+ BioRestartLong = 'DA__FREQ_1/RST_after.'//DA_DATE//'.'//DA_VarList(iVar)//'.nc'
+
+ if(drv%Verbose .eq. 1 .and. MyId .eq. 0) &
+ print*, "Writing Nut Restart ", BioRestart
+
+ ierr = nf90mpi_create(Var3DCommunicator, BioRestart, NF90_CLOBBER, MPI_INFO_NULL, ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_create '//BioRestart, ierr)
+
+ ierr = nf90mpi_def_dim(ncid,'x',global_im ,xid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim longitude ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'y' ,global_jm ,yid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim latitude ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'z' ,global_km, depid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim depth ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'time',MyTime ,timeId)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim time ', ierr)
+
+ MyVarName='TRN'//DA_VarList(iVar)
+
+ ierr = nf90mpi_def_var(ncid, MyVarName, nf90_float, (/xid,yid,depid,timeId/), idP )
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', ierr)
+
+ ierr = nf90mpi_def_var(ncid,'time' , nf90_double, (/timeid/) , idTim)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', ierr)
+ ierr = nf90mpi_put_att(ncid,idP , 'missing_value',1.e+20)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_att', ierr)
+
+ ierr = nf90mpi_enddef(ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_enddef'//DA_VarList(iVar), ierr)
+
+ do k=1,grd%km
+ do j=1,grd%jm
+ do i=1,grd%im
+
+ if(bio%InitialNut(i,j,k,1) .lt. 1.e20) then
+ DumpBio(i,j,k,1) = bio%InitialNut(i,j,k,l)
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ ierr = nf90mpi_put_var_all(ncid,idP,DumpBio,MyStart_4d,MyCount_4d)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all '//DA_VarList(iVar), ierr)
+
+ ierr = nf90mpi_put_var_all(ncid,idTim,TimeArr,MyStartSingle,MyCountSingle)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all '//DA_VarList(iVar), ierr)
+
+ ierr = nf90mpi_close(ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_close '//BioRestart, ierr)
+
+ call MPI_Barrier(Var3DCommunicator, ierr)
+ ! only process 0 creates link to restart files
+ if(MyId .eq. 0) then
+ SysErr = system("ln -sf $PWD/"//BioRestart//" "//BioRestartLong)
+ if(SysErr /= 0) call MPI_Abort(MPI_COMM_WORLD, -1, SysErr)
+ endif
+ enddo ! l
+
+ DEALLOCATE(DumpBio)
+ ! DEALLOCATE(DumpBio, ValuesToTest, MyConditions)
+
+end subroutine cp_nut_stat
diff --git a/da_params.f90 b/da_params.f90
index 9cc98e8..28f6622 100644
--- a/da_params.f90
+++ b/da_params.f90
@@ -11,16 +11,16 @@ MODULE DA_PARAMS
integer :: NNutVar ! number of nutrient variables
integer :: NBioVar ! number of nutrient variables
CHARACTER(LEN=3), allocatable :: DA_VarList(:) ! name of DA biological variables
- integer :: DA_JulianDate ! julian date
+ double precision :: DA_JulianDate ! julian date
CONTAINS
SUBROUTINE SET_DA_PARAMS
- DA_DATE = '20150101-12:00:00'
+ DA_DATE = '20150101-00:00:00'
ShortDate = DA_DATE(1:11)//DA_DATE(13:14)//DA_DATE(16:17)
- jpk_200 = 36
- NPhytoVar = 0
+ jpk_200 = 60
+ NPhytoVar = 17
NNutVar = 2
NBioVar = NPhytoVar + NNutVar
@@ -29,33 +29,33 @@ SUBROUTINE SET_DA_PARAMS
! DA_VarList init
! It must be consistent with NPhytoVar and NNutVar values
- ! DA_VarList( 1)='P1l'
- ! DA_VarList( 2)='P2l'
- ! DA_VarList( 3)='P3l'
- ! DA_VarList( 4)='P4l'
+ DA_VarList( 1)='P1l'
+ DA_VarList( 2)='P2l'
+ DA_VarList( 3)='P3l'
+ DA_VarList( 4)='P4l'
- ! DA_VarList( 5)='P1c'
- ! DA_VarList( 6)='P2c'
- ! DA_VarList( 7)='P3c'
- ! DA_VarList( 8)='P4c'
+ DA_VarList( 5)='P1c'
+ DA_VarList( 6)='P2c'
+ DA_VarList( 7)='P3c'
+ DA_VarList( 8)='P4c'
- ! DA_VarList( 9)='P1n'
- ! DA_VarList(10)='P2n'
- ! DA_VarList(11)='P3n'
- ! DA_VarList(12)='P4n'
+ DA_VarList( 9)='P1n'
+ DA_VarList(10)='P2n'
+ DA_VarList(11)='P3n'
+ DA_VarList(12)='P4n'
- ! DA_VarList(13)='P1p'
- ! DA_VarList(14)='P2p'
- ! DA_VarList(15)='P3p'
- ! DA_VarList(16)='P4p'
+ DA_VarList(13)='P1p'
+ DA_VarList(14)='P2p'
+ DA_VarList(15)='P3p'
+ DA_VarList(16)='P4p'
- ! DA_VarList(17)='P1s'
+ DA_VarList(17)='P1s'
- ! DA_VarList(18)='N3n'
- ! DA_VarList(19)='N1p'
+ DA_VarList(18)='N3n'
+ DA_VarList(19)='N1p'
- DA_VarList(1)='N3n'
- DA_VarList(2)='N1p'
+ ! DA_VarList(1)='N3n'
+ ! DA_VarList(2)='N1p'
END SUBROUTINE SET_DA_PARAMS
@@ -65,4 +65,4 @@ SUBROUTINE CLEAN_DA_PARAMS
END SUBROUTINE CLEAN_DA_PARAMS
-END MODULE DA_PARAMS
\ No newline at end of file
+END MODULE DA_PARAMS
diff --git a/def_cov.f90 b/def_cov.f90
index a775e34..de48ad7 100644
--- a/def_cov.f90
+++ b/def_cov.f90
@@ -36,10 +36,12 @@ subroutine def_cov
use rcfl
use mpi_str
use bio_str
+ use da_params
implicit none
INTEGER(i4) :: k, nspl, i, j, kk, l
+ INTEGER(i4) :: nsplvec(1)
REAL(r8) :: E, dst, Lmean, mean_rad
REAL(r8) , ALLOCATABLE :: sfct(:), al(:), bt(:)
INTEGER(i4) , ALLOCATABLE :: jnxx(:)
@@ -108,6 +110,7 @@ subroutine def_cov
!nspl = max(grd%jm,grd%im)
nspl = max(GlobalRow,GlobalCol)
+ nsplvec(:) = nspl
ALLOCATE ( sfct(nspl)) ; sfct = huge(sfct(1))
ALLOCATE ( jnxx(nspl)) ; jnxx = huge(jnxx(1))
ALLOCATE ( al(nspl)) ; al = huge(al(1))
@@ -154,8 +157,8 @@ subroutine def_cov
jnxx(j) = j
enddo
sfct(nspl/2+1) = 1.
- call rcfl_y_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nspl)
- call rcfl_y_ad_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nspl)
+ call rcfl_y_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nsplvec)
+ call rcfl_y_ad_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nsplvec)
rcf%sc(k,l) = sfct(nspl/2+1)
enddo
enddo
@@ -232,8 +235,8 @@ subroutine def_cov
jnxx(j) = j
enddo
sfct(nspl/2+1) = 1.
- call rcfl_y_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nspl)
- call rcfl_y_ad_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nspl)
+ call rcfl_y_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nsplvec)
+ call rcfl_y_ad_init ( 1, nspl, 1, nspl, al, bt, sfct, jnxx, nsplvec)
rcf%sc(k,l) = sfct(nspl/2+1)
enddo
enddo
@@ -330,6 +333,57 @@ subroutine def_cov
enddo
+
+
+ if (drv%multiv.eq.1) then
+ grd%imax_chl = 0
+ grd%jmax_chl = 0
+
+ do k = 1, ros%kmchl
+
+ grd%imx(k) = 0
+ do j = 1, localCol
+ kk = grd%istp(1,j,k)
+ if( grd%global_msk(1,j+GlobalColOffset,k).eq.1. ) kk = kk + 1
+ grd%inx(1,j,k) = kk
+ do i = 2, GlobalRow
+ if( grd%global_msk(i,j+GlobalColOffset,k).eq.0. .and. grd%global_msk(i-1,j+GlobalColOffset,k).eq.1. ) then
+ kk = kk + grd%istp(i,j,k)
+ else if( grd%global_msk(i,j+GlobalColOffset,k).eq.1. .and. grd%global_msk(i-1,j+GlobalColOffset,k).eq.0. ) then
+ kk = kk + grd%istp(i,j,k) + 1
+ else if( grd%global_msk(i,j+GlobalColOffset,k).eq.1. ) then
+ kk = kk + 1
+ endif
+ grd%inx(i,j,k) = kk
+ enddo
+ grd%imx(k) = max( grd%imx(k), kk+grd%istp(GlobalRow,j,k))
+ enddo
+ grd%imax = max( grd%imax, grd%imx(k))
+
+ grd%jmx(k) = 0
+ do i = 1, localRow
+ kk = grd%jstp(i,1,k)
+ if( grd%global_msk(i+GlobalRowOffset,1,k).eq.1. ) kk = kk + 1
+ grd%jnx(i,1,k) = kk
+ do j = 2, GlobalCol
+ if( grd%global_msk(i+GlobalRowOffset,j,k).eq.0. .and. grd%global_msk(i+GlobalRowOffset,j-1,k).eq.1. ) then
+ kk = kk + grd%jstp(i,j,k)
+ else if( grd%global_msk(i+GlobalRowOffset,j,k).eq.1. .and. grd%global_msk(i+GlobalRowOffset,j-1,k).eq.0. ) then
+ kk = kk + grd%jstp(i,j,k) + 1
+ else if( grd%global_msk(i+GlobalRowOffset,j,k).eq.1. ) then
+ kk = kk + 1
+ endif
+ grd%jnx(i,j,k) = kk
+ enddo
+ grd%jmx(k) = max( grd%jmx(k), kk+grd%jstp(i,GlobalCol,k))
+ enddo
+ grd%jmax = max( grd%jmax, grd%jmx(k))
+
+ enddo
+ call MPI_Allreduce(MPI_IN_PLACE, grd%imax_chl, 1, MPI_INT, MPI_MAX, Var3DCommunicator, ierr)
+ call MPI_Allreduce(MPI_IN_PLACE, grd%jmax_chl, 1, MPI_INT, MPI_MAX, Var3DCommunicator, ierr)
+ endif ! if drv%multiv
+
call MPI_Allreduce(MPI_IN_PLACE, grd%imax, 1, MPI_INT, MPI_MAX, Var3DCommunicator, ierr)
call MPI_Allreduce(MPI_IN_PLACE, grd%jmax, 1, MPI_INT, MPI_MAX, Var3DCommunicator, ierr)
@@ -396,28 +450,40 @@ subroutine def_cov
ros%kmt = grd%km
- if(drv%chl_assim .eq. 1) then
- call rdeofs_chl
- else
- ros%neof_chl = 0
- endif
+ if (drv%multiv .eq. 0) then
+ ros%neof_multi = 0
- if(drv%nut .eq. 1) then
- if(bio%n3n .eq. 1) then
- call rdeofs_n3n
+ if(drv%chl_assim .eq. 1) then
+ call rdeofs_chl
else
- ros%neof_n3n = 0
+ ros%neof_chl = 0
endif
- if(bio%o2o .eq. 1) then
- call rdeofs_o2o
+
+ if(drv%nut .eq. 1) then
+ if(bio%n3n .eq. 1) then
+ call rdeofs_n3n
+ else
+ ros%neof_n3n = 0
+ endif
+ if(bio%o2o .eq. 1) then
+ call rdeofs_o2o
+ else
+ ros%neof_o2o = 0
+ endif
+ ros%neof_nut = ros%neof_n3n + ros%neof_o2o
else
- ros%neof_o2o = 0
+ ros%neof_nut = 0
endif
- ros%neof_nut = ros%neof_n3n + ros%neof_o2o
- else
- ros%neof_nut = 0
+
+ else if(drv%multiv .eq. 1) then
+ ros%neof_chl = 0
+ ros%neof_n3n = 0
+ ros%neof_o2o = 0
+ call rdeofs_multi
endif
- ros%neof = ros%neof_chl + ros%neof_nut
+
+
+ ros%neof = ros%neof_multi + ros%neof_chl + ros%neof_nut
ALLOCATE ( grd%ro( grd%im, grd%jm, ros%neof)) ; grd%ro = 0.0
ALLOCATE ( grd%ro_ad( grd%im, grd%jm, ros%neof)) ; grd%ro_ad = 0.0
@@ -425,6 +491,11 @@ subroutine def_cov
if(MyId .eq. 0) then
write(*,*) 'rcfl allocation :', grd%jmax, grd%imax, nthreads
write(*,*) 'Total number of eofs: ', ros%neof
+ write(*,*) ' - multi number of eofs: ', ros%neof_multi
+ write(*,*) ' - chl number of eofs: ', ros%neof_chl
+ write(*,*) ' - nut number of eofs: ', ros%neof_nut
+ write(*,*) ' - n3n number of eofs: ', ros%neof_n3n
+ write(*,*) ' - o2o number of eofs: ', ros%neof_o2o
endif
ALLOCATE ( a_rcx(localCol,grd%imax,nthreads)) ; a_rcx = huge(a_rcx(1,1,1))
ALLOCATE ( b_rcx(localCol,grd%imax,nthreads)) ; b_rcx = huge(b_rcx(1,1,1))
@@ -447,15 +518,30 @@ subroutine def_cov
! read ratios for biological repartition
! of the chlorophyll
- if(drv%chl_assim.eq.1) then
- call readChlStat
- endif
+ if(drv%multiv .eq. 0) then
+ if(drv%chl_assim.eq.1) then
+ call readChlStat
+ if ((drv%nut .eq. 0) .and. (NNutVar .gt. 0)) then
+ call readNutStat
+ if (drv%chl_upnut.eq.1) then
+ call readNutCov
+ call readChlNutCov
+ endif
+ endif
+ endif
- if(drv%nut.eq.1) then
- call readNutStat
- if(bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) then
- call readNutCov
- endif
- endif
+ if(drv%nut.eq.1) then
+ call readNutStat
+ if(bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) then
+ call readNutCov
+ endif
+ endif
+
+ else if (drv%multiv.eq.1) then
+ call readChlStat
+ call readNutStat
+ if(bio%updateN1p.eq.1) &
+ call readNutCov
+ endif
end subroutine def_cov
diff --git a/def_nml.f90 b/def_nml.f90
index c0cbf33..f99c065 100644
--- a/def_nml.f90
+++ b/def_nml.f90
@@ -40,17 +40,14 @@ subroutine def_nml
implicit none
- LOGICAL :: read_eof, ApplyConditions
+ LOGICAL :: read_eof
INTEGER(i4) :: neof_chl, neof_n3n, neof_o2o, nreg, rcf_ntr
- INTEGER(i4) :: ctl_m, chl_assim, nut, N3n, O2o, updateN1p
- INTEGER(i4) :: biol, bphy, nphyto, uniformL, anisL, verbose
- REAL(r8) :: rcf_L, ctl_tol, ctl_per, rcf_efc, chl_dep
- INTEGER(i4) :: argo, sat_obs, ncmp
+ INTEGER(i4) :: neof_multi, kmchl, kmnit
+ INTEGER(i4) :: verbose
+ REAL(r8) :: rcf_L, ctl_tol, ctl_per, rcf_efc
- NAMELIST /ctllst/ ctl_tol, ctl_per
- NAMELIST /covlst/ neof_chl, neof_n3n, neof_o2o, nreg, read_eof, rcf_ntr, rcf_L, rcf_efc
- NAMELIST /biolst/ chl_assim, nut, nphyto, chl_dep, ncmp, ApplyConditions, N3n, updateN1p, O2o
- NAMELIST /params/ sat_obs, argo, uniformL, anisL, verbose
+ NAMELIST /ctllst/ ctl_tol, ctl_per, verbose
+ NAMELIST /covlst/ neof_chl, neof_n3n, neof_o2o, neof_multi, kmchl, kmnit, nreg, read_eof, rcf_ntr, rcf_L, rcf_efc
! -------------------------------------------------------------------
@@ -60,7 +57,7 @@ subroutine def_nml
drv%dia = 12
if(MyId .eq. 0) then
- open ( drv%dia, file='OceanVar.diagnostics', form='formatted' )
+ open ( drv%dia, file='BioVar.diagnostics', form='formatted' )
endif
!---------------------------------------------------------------------
@@ -83,11 +80,13 @@ subroutine def_nml
write(drv%dia,*) ' MINIMIZER NAMELIST INPUT: '
write(drv%dia,*) ' Minimum gradient of J: ctl_tol = ', ctl_tol
write(drv%dia,*) ' Percentage of initial gradient: ctl_per = ', ctl_per
+ write(drv%dia,*) ' Add verbose on standard output verbose = ', verbose
endif
ctl%pgtol = ctl_tol
ctl%pgper = ctl_per
+ drv%Verbose = verbose
! ---
read(11,covlst)
@@ -100,6 +99,9 @@ subroutine def_nml
write(drv%dia,*) ' Number of EOFs for chl: neof_chl = ', neof_chl
write(drv%dia,*) ' Number of EOFs for N3n: neof_n3n = ', neof_n3n
write(drv%dia,*) ' Number of EOFs for O2o: neof_o2o = ', neof_o2o
+ write(drv%dia,*) ' Number of multivariate EOFs: neof_multi = ', neof_multi
+ write(drv%dia,*) ' Chl Nlevels in multi EOFs: kmchl = ', kmchl
+ write(drv%dia,*) ' Nit Nlevels in multi EOFs: kmnit = ', kmnit
write(drv%dia,*) ' Number of regions: nreg = ', nreg
write(drv%dia,*) ' Read EOFs from a file: read_eof = ', read_eof
write(drv%dia,*) ' Half number of iterations: rcf_ntr = ', rcf_ntr
@@ -108,70 +110,19 @@ subroutine def_nml
endif
+ close(11)
+
ros%neof_chl = neof_chl
ros%neof_n3n = neof_n3n
ros%neof_o2o = neof_o2o
+ ros%neof_multi = neof_multi
+ ros%kmchl = kmchl
+ ros%kmnit = kmnit
ros%nreg = nreg
ros%read_eof = read_eof
rcf%ntr = rcf_ntr
rcf%L = rcf_L
rcf%efc = rcf_efc
-! ---
- read(11,biolst)
-
- if(MyId .eq. 0) then
-
- write(drv%dia,*) '------------------------------------------------------------'
- write(drv%dia,*) '------------------------------------------------------------'
- write(drv%dia,*) ' BIOLOGY NAMELIST INPUT: '
- write(drv%dia,*) ' Chlorophyll assimilation chl_assim = ', chl_assim
- write(drv%dia,*) ' Nutrient assimilation nut = ', nut
- write(drv%dia,*) ' Number of phytoplankton species nphyt = ', nphyto
- write(drv%dia,*) ' Minimum depth for chlorophyll chl_dep = ', chl_dep
- write(drv%dia,*) ' Number of phytoplankton components ncmp = ', ncmp
- write(drv%dia,*) ' Apply conditions flag ApplyConditions = ', ApplyConditions
- write(drv%dia,*) ' N3n assimilation N3n = ', N3n
- write(drv%dia,*) ' N1p update based on N3n assimilation updateN1p = ', updateN1p
- write(drv%dia,*) ' O2o assimilation O2o = ', O2o
-
- endif
-
- drv%chl_assim = chl_assim
- drv%nut = nut
- bio%nphy = nphyto
- sat%dep = chl_dep
- bio%ncmp = ncmp
- bio%ApplyConditions = ApplyConditions
- bio%N3n = N3n
- bio%updateN1p = updateN1p
- bio%O2o = O2o
-
- read(11,params)
-
- if(MyId .eq. 0) then
-
- write(drv%dia,*) '------------------------------------------------------------'
- write(drv%dia,*) '------------------------------------------------------------'
- write(drv%dia,*) ' PARAMETERS NAMELIST INPUT: '
- write(drv%dia,*) ' Read Satellite observations sat_obs = ', sat_obs
- write(drv%dia,*) ' Read ARGO float observations argo = ', argo
- write(drv%dia,*) ' Set uniform correlation radius uniformL = ', uniformL
- write(drv%dia,*) ' Set anisotropy on corr radius anisL = ', anisL
- write(drv%dia,*) ' Add verbose on standard output verbose = ', verbose
- write(drv%dia,*) '------------------------------------------------------------'
- write(drv%dia,*) '------------------------------------------------------------'
- write(drv%dia,*) ''
-
-
- endif
-
- close(11)
-
- drv%sat_obs = sat_obs
- drv%argo_obs = argo
- drv%uniformL = uniformL
- drv%anisL = anisL
- drv%Verbose = verbose
end subroutine def_nml
diff --git a/def_nml_multi.f90 b/def_nml_multi.f90
new file mode 100644
index 0000000..d5a0293
--- /dev/null
+++ b/def_nml_multi.f90
@@ -0,0 +1,138 @@
+subroutine def_nml_multi
+
+!---------------------------------------------------------------------------
+! !
+! Copyright 2018 Anna Teruzzi, OGS, Trieste !
+! !
+! This file is part of 3DVarBio. !
+! 3DVarBio is based on OceanVar (Dobricic, 2006) !
+! !
+! 3DVarBio is free software: you can redistribute it and/or modify. !
+! it under the terms of the GNU General Public License as published by !
+! the Free Software Foundation, either version 3 of the License, or !
+! (at your option) any later version. !
+! !
+! 3DVarBio is distributed in the hope that it will be useful, !
+! but WITHOUT ANY WARRANTY; without even the implied warranty of !
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
+! GNU General Public License for more details. !
+! !
+! You should have received a copy of the GNU General Public License !
+! along with OceanVar. If not, see . !
+! !
+!---------------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! !
+! Define analysis parameters from namelists !
+! for multi platform and multivariate DA !
+! (other general DA parameters are defined in def_nml.f90) !
+! !
+! Version 1: A. Teruzzi 2018 !
+!-----------------------------------------------------------------------
+
+
+ use set_knd
+ use drv_str
+ use grd_str
+ use obs_str
+ use eof_str
+ use cns_str
+ use ctl_str
+ use mpi_str
+ use bio_str
+ use da_params
+
+ implicit none
+
+ LOGICAL :: ApplyConditions
+ INTEGER(i4) :: chl_assim, chl_upnut, nut, multiv, N3n, O2o, updateN1p
+ INTEGER(i4) :: nphyto, uniformL, anisL
+ REAL(r8) :: chl_dep
+ INTEGER(i4) :: argo, sat_obs, ncmp
+
+ !NAMELIST /ctllst/ ctl_tol, ctl_per
+ !NAMELIST /covlst/ neof_chl, neof_n3n, neof_o2o, nreg, read_eof, rcf_ntr, rcf_L, rcf_efc
+ NAMELIST /biolst/ chl_assim, chl_upnut, nut, multiv, nphyto, chl_dep, ncmp, ApplyConditions, N3n, updateN1p, O2o
+ NAMELIST /params/ sat_obs, argo, uniformL, anisL
+
+
+! -------------------------------------------------------------------
+! Open a formatted file for the diagnostics
+! ---
+
+ drv%dia = 12
+
+ if(MyId .eq. 0) then
+ open ( drv%dia, file='DA__FREQ_1/OceanVar.dia_multinml.'//DA_DATE, form='formatted' )
+ endif
+
+!---------------------------------------------------------------------
+! Open the namelist
+! ---
+
+ open(11,file='DA__FREQ_1/satfloat.'//DA_DATE//'.nml',form='formatted')
+
+ ! ---
+
+ read(11,biolst)
+
+ if(MyId .eq. 0) then
+
+ write(drv%dia,*) '------------------------------------------------------------'
+ write(drv%dia,*) '------------------------------------------------------------'
+ write(drv%dia,*) ' BIOLOGY NAMELIST INPUT: '
+ write(drv%dia,*) ' Chlorophyll assimilation chl_assim = ', chl_assim
+ write(drv%dia,*) ' N3n update based on chl assimilation chl_upnut = ', chl_upnut
+ write(drv%dia,*) ' Nutrient assimilation nut = ', nut
+ write(drv%dia,*) ' Multivariate assimilation multiv = ', multiv
+ write(drv%dia,*) ' Number of phytoplankton species nphyt = ', nphyto
+ write(drv%dia,*) ' Minimum depth for chlorophyll chl_dep = ', chl_dep
+ write(drv%dia,*) ' Number of phytoplankton components ncmp = ', ncmp
+ write(drv%dia,*) ' Apply conditions flag ApplyConditions = ', ApplyConditions
+ write(drv%dia,*) ' N3n assimilation N3n = ', N3n
+ write(drv%dia,*) ' N1p update based on N3n assimilation updateN1p = ', updateN1p
+ write(drv%dia,*) ' O2o assimilation O2o = ', O2o
+
+ endif
+
+ drv%chl_assim = chl_assim
+ drv%chl_upnut = chl_upnut
+ drv%nut = nut
+ drv%multiv = multiv
+ bio%nphy = nphyto
+ sat%dep = chl_dep
+ bio%ncmp = ncmp
+ bio%ApplyConditions = ApplyConditions
+ bio%N3n = N3n
+ bio%updateN1p = updateN1p
+ bio%O2o = O2o
+
+ ! ---
+
+ read(11,params)
+
+ if(MyId .eq. 0) then
+
+ write(drv%dia,*) '------------------------------------------------------------'
+ write(drv%dia,*) '------------------------------------------------------------'
+ write(drv%dia,*) ' PARAMETERS NAMELIST INPUT: '
+ write(drv%dia,*) ' Read Satellite observations sat_obs = ', sat_obs
+ write(drv%dia,*) ' Read ARGO float observations argo = ', argo
+ write(drv%dia,*) ' Set uniform correlation radius uniformL = ', uniformL
+ write(drv%dia,*) ' Set anisotropy on corr radius anisL = ', anisL
+ write(drv%dia,*) '------------------------------------------------------------'
+ write(drv%dia,*) '------------------------------------------------------------'
+ write(drv%dia,*) ''
+
+
+ endif
+
+ close(11)
+
+ drv%sat_obs = sat_obs
+ drv%argo_obs = argo
+ drv%uniformL = uniformL
+ drv%anisL = anisL
+
+end subroutine def_nml_multi
diff --git a/drv_str.f90 b/drv_str.f90
index 9c3e616..3e73d84 100644
--- a/drv_str.f90
+++ b/drv_str.f90
@@ -42,10 +42,12 @@ MODULE drv_str
INTEGER(i4) :: sat_obs ! Flag for the assimilation of the satellite observations
INTEGER(i4) :: argo_obs ! Flag for the assimilation of the argo float observations
INTEGER(i4) :: chl_assim ! Flag for the chlorophyll assimilation
+ INTEGER(i4) :: chl_upnut ! Flag for the update of nut based on chlorophyll assimilation
INTEGER(i4) :: uniformL ! Flag for setting uniform correlation radius (1 = non uniform)
INTEGER(i4) :: anisL ! Flag for setting anisotropy on correlation radius (1 = anisotropy)
INTEGER(i4) :: Verbose ! Flag for printing verbose output
- INTEGER(i4) :: nut ! Flag for the chlorophyll assimilation
+ INTEGER(i4) :: nut ! Flag for nutrient assimilation
+ INTEGER(i4) :: multiv ! Flag for multivariate assimilation
END TYPE drv_t
diff --git a/eof_str.f90 b/eof_str.f90
index bbbba9f..75cd61c 100644
--- a/eof_str.f90
+++ b/eof_str.f90
@@ -42,8 +42,11 @@ MODULE eof_str
INTEGER(i4) :: neof_nut ! No. of EOFs for nutrients
INTEGER(i4) :: neof_n3n ! No. of EOFs for N3n
INTEGER(i4) :: neof_o2o ! No. of EOFs for O2o
+ INTEGER(i4) :: neof_multi ! No. of EOFs for multivariate
INTEGER(i4) :: nreg ! No. of regions
INTEGER(i4) :: kmt ! No. of levels of EOFs
+ INTEGER(i4) :: kmchl ! No. of levels of multi EOFs for chl
+ INTEGER(i4) :: kmnit ! No. of levels of multi EOFs for nit
REAL(r8), POINTER :: evcr(:,:,:) ! Eigenvectors on regions
REAL(r8), POINTER :: evar(:,:) ! Eigenvalues on regions
REAL(r8), POINTER :: corr(:,:,:) ! Corelations on regions
@@ -57,6 +60,8 @@ MODULE eof_str
REAL(r8), POINTER :: eva_n3n(:,:) ! Eigenvalues
REAL(r8), POINTER :: evc_o2o(:,:,:) ! Eigenvectors
REAL(r8), POINTER :: eva_o2o(:,:) ! Eigenvalues
+ REAL(r8), POINTER :: evc_multi(:,:,:) ! Eigenvectors
+ REAL(r8), POINTER :: eva_multi(:,:) ! Eigenvalues
#endif
diff --git a/filename_mod.f90 b/filename_mod.f90
index 8eb8c8d..dcf6e14 100644
--- a/filename_mod.f90
+++ b/filename_mod.f90
@@ -5,8 +5,11 @@ MODULE FILENAMES
character (LEN=1024) :: EOF_FILE_CHL != 'eofs_chl.nc'
character (LEN=1024) :: EOF_FILE_N3N != 'eofs_n3n.nc'
character (LEN=1024) :: EOF_FILE_O2O != 'eofs_o2o.nc'
+character (LEN=1024) :: EOF_FILE_MULTI != 'eofs_multi.nc'
+character (LEN=1024) :: STD_FILE_MULTI != 'std_multi.nc'
character (LEN=1024) :: MISFIT_FILE != 'chl_mis.nc'
-character (LEN=1024) :: NUTCOV_FILE != 'cov_N3nN1p.nc'
+character (LEN=1024) :: NUTCOV_FILE != 'crosscorrs.nc'
+character (LEN=1024) :: NUTCHLCOV_FILE != 'crosscorrs.nc'
character (LEN=1024) :: CORR_FILE != 'corr.nc'
character (LEN=1024) :: EIV_FILE != 'eiv.nc'
character (LEN=1024) :: OBS_FILE != 'obs_1.dat'
@@ -21,13 +24,16 @@ MODULE FILENAMES
!
!
SUBROUTINE SETFILENAMES
-! !VAR_FILE='DA_static_data/MISFIT/VAR2D/var2D.05.nc'
+! !VAR_FILE='DA_static_data/VAR_SAT/var2D.05.nc'
!
EOF_FILE_CHL = 'eofs_chl.nc'
EOF_FILE_N3N = 'eofs_n3n.nc'
EOF_FILE_O2O = 'eofs_o2o.nc'
+EOF_FILE_MULTI = 'eofs_multi.nc'
+STD_FILE_MULTI = 'std_multi.nc'
MISFIT_FILE = 'chl_mis.nc'
-NUTCOV_FILE = 'cov_N3nN1p.nc'
+NUTCOV_FILE = 'crosscorrs.nc'
+NUTCHLCOV_FILE = 'crosscorrs.nc'
CORR_FILE = 'corr.nc'
EIV_FILE = 'eiv.nc'
OBS_FILE = 'obs_1.dat' ! 'obs_'//fgrd//'.dat'
diff --git a/get_obs.f90 b/get_obs.f90
index c87cbc0..ed087a3 100644
--- a/get_obs.f90
+++ b/get_obs.f90
@@ -37,6 +37,8 @@ subroutine get_obs
arg%no = 0
sat%no = 0
+ arg%nc = 0
+ sat%nc = 0
! ----
diff --git a/get_obs_arg.f90 b/get_obs_arg.f90
index 24da036..7a45d6d 100644
--- a/get_obs_arg.f90
+++ b/get_obs_arg.f90
@@ -41,6 +41,7 @@ subroutine get_obs_arg
INTEGER(i4) :: k
INTEGER(i4) :: i1, kk, i
REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpFlc, TmpPar, TmpLon, TmpLat
+ ! REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpDpt, TmpTim, TmpRes, TmpErr, TmpStd, TmpIno
REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpDpt, TmpTim, TmpRes, TmpErr, TmpIno
INTEGER(i4) :: GlobalArgNum, Counter, ierr
character(len=1024) :: filename
@@ -69,6 +70,7 @@ subroutine get_obs_arg
ALLOCATE( TmpLon(GlobalArgNum), TmpLat(GlobalArgNum))
ALLOCATE( TmpDpt(GlobalArgNum), TmpTim(GlobalArgNum))
ALLOCATE( TmpRes(GlobalArgNum), TmpErr(GlobalArgNum))
+! ALLOCATE( TmpStd(GlobalArgNum), TmpIno(GlobalArgNum))
ALLOCATE( TmpIno(GlobalArgNum))
if(MyId .eq. 0) then
@@ -78,7 +80,9 @@ subroutine get_obs_arg
TmpFlc(k), TmpPar(k), &
TmpLon(k), TmpLat(k), &
TmpDpt(k), TmpTim(k), &
- TmpRes(k), TmpErr(k), TmpIno(k)
+ TmpRes(k), TmpErr(k), &
+ ! TmpStd(k), TmpIno(k)
+ TmpIno(k)
end do
close (511)
endif
@@ -98,9 +102,9 @@ subroutine get_obs_arg
do k=1,GlobalArgNum
if( TmpLon(k) .ge. grd%lon(1,1) .and. TmpLon(k) .lt. grd%NextLongitude .and. &
TmpLat(k) .ge. grd%lat(1,1) .and. TmpLat(k) .lt. grd%lat(grd%im,grd%jm) ) then
- if((TmpPar(k).eq.0 .and. drv%chl_assim.eq.1) .or. &
- (TmpPar(k).eq.1 .and. drv%nut.eq.1 .and. bio%n3n.eq.1) .or. &
- (TmpPar(k).eq.2 .and. drv%nut.eq.1 .and. bio%o2o.eq.1)) then
+ if( (TmpPar(k).eq.0 .and. ((drv%chl_assim.eq.1).or.(drv%multiv.eq.1))) .or. &
+ (TmpPar(k).eq.1 .and. ((drv%nut.eq.1 .and.bio%n3n.eq.1).or.(drv%multiv.eq.1))) .or. &
+ (TmpPar(k).eq.2 .and. drv%nut.eq.1 .and. bio%o2o.eq.1) ) then
Counter = Counter + 1
endif
endif
@@ -116,6 +120,7 @@ subroutine get_obs_arg
ALLOCATE ( arg%inc(arg%no))
ALLOCATE ( arg%err(arg%no))
ALLOCATE ( arg%res(arg%no))
+! ALLOCATE ( arg%std(arg%no))
ALLOCATE ( arg%ib(arg%no), arg%jb(arg%no), arg%kb(arg%no))
ALLOCATE ( arg%pb(arg%no), arg%qb(arg%no), arg%rb(arg%no))
ALLOCATE ( arg%pq1(arg%no), arg%pq2(arg%no), arg%pq3(arg%no), arg%pq4(arg%no))
@@ -125,8 +130,8 @@ subroutine get_obs_arg
do k=1,GlobalArgNum
if( TmpLon(k) .ge. grd%lon(1,1) .and. TmpLon(k) .lt. grd%NextLongitude .and. &
TmpLat(k) .ge. grd%lat(1,1) .and. TmpLat(k) .lt. grd%lat(grd%im,grd%jm) ) then
- if((TmpPar(k).eq.0 .and. drv%chl_assim.eq.1) .or. &
- (TmpPar(k).eq.1 .and. drv%nut.eq.1 .and. bio%n3n.eq.1) .or. &
+ if((TmpPar(k).eq.0 .and. (drv%chl_assim.eq.1 .or. drv%multiv.eq.1)) .or. &
+ (TmpPar(k).eq.1 .and. ((drv%nut.eq.1 .and. bio%n3n.eq.1).or.(drv%multiv.eq.1))) .or. &
(TmpPar(k).eq.2 .and. drv%nut.eq.1 .and. bio%o2o.eq.1)) then
Counter = Counter + 1
arg%flc(Counter) = TmpFlc(k)
@@ -136,6 +141,7 @@ subroutine get_obs_arg
arg%dpt(Counter) = TmpDpt(k)
arg%res(Counter) = TmpRes(k)
arg%err(Counter) = TmpErr(k)
+ ! arg%std(Counter) = TmpStd(k)
arg%ino(Counter) = TmpIno(k)
endif
endif
@@ -194,6 +200,7 @@ subroutine get_obs_arg
DEALLOCATE( TmpLon, TmpLat)
DEALLOCATE( TmpDpt, TmpTim)
DEALLOCATE( TmpRes, TmpErr)
+! DEALLOCATE( TMPStd, TmpIno)
DEALLOCATE( TmpIno)
end subroutine get_obs_arg
@@ -212,12 +219,13 @@ subroutine int_par_arg
use set_knd
use drv_str
use grd_str
+ use eof_str
use obs_str
use mpi_str
implicit none
- INTEGER(i4) :: i, j, k, ierr
+ INTEGER(i4) :: i, j, k, ierr, kind, kk
INTEGER(i4) :: i1, j1, k1, idep
REAL(r8) :: p1, q1, r1
REAL(r8) :: msk4, div_x, div_y
@@ -280,7 +288,8 @@ subroutine int_par_arg
! ---
! Horizontal interpolation parameters for each masked grid
do k = 1,arg%no
- if(arg%flc(k) .eq. 1) then
+ if(arg%flg(k) .eq. 1) then
+ ! if(arg%flg(k) .eq. 1) then ! to verify that it works also in this case
i1=arg%ib(k)
p1=arg%pb(k)
@@ -355,6 +364,25 @@ subroutine int_par_arg
endif
enddo
+
+ ! Exclude observations below ros%kmchl in multivariate observations
+ if(drv%multiv.eq.1) then
+ do k=1,arg%no
+ if((arg%flc(k).eq.1).and.(arg%par(k).eq.0)) then
+ kind = grd%km-1
+ do kk = 1,grd%km-1
+ if( arg%dpt(k).ge.grd%dep(kk) .and. arg%dpt(k).lt.grd%dep(kk+1) ) then
+ kind = kk
+ else if ( arg%dpt(k).ge.0 .and. arg%dpt(k).lt.grd%dep(1)) then
+ kind = 1
+ endif
+ enddo
+ if(kind.gt.ros%kmchl) then
+ arg%flc(k)=0
+ end if
+ endif
+ enddo
+ endif
! ---
! Count good observations
diff --git a/get_obs_sat.f90 b/get_obs_sat.f90
index 250730d..0cc426f 100644
--- a/get_obs_sat.f90
+++ b/get_obs_sat.f90
@@ -2,16 +2,17 @@ subroutine get_obs_sat
!---------------------------------------------------------------------------
! !
- ! Copyright 2006 Srdjan Dobricic, CMCC, Bologna !
+ ! Copyright 2018 Anna Teruzzi, OGS, Trieste !
! !
- ! This file is part of OceanVar. !
+ ! This file is part of 3DVarBio.
+ ! 3DVarBio is based on OceanVar (Dobricic, 2006) !
! !
- ! OceanVar is free software: you can redistribute it and/or modify. !
+ ! 3DVarBio is free software: you can redistribute it and/or modify. !
! it under the terms of the GNU General Public License as published by !
! the Free Software Foundation, either version 3 of the License, or !
! (at your option) any later version. !
! !
- ! OceanVar is distributed in the hope that it will be useful, !
+ ! 3DVarBio is distributed in the hope that it will be useful, !
! but WITHOUT ANY WARRANTY; without even the implied warranty of !
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
! GNU General Public License for more details. !
@@ -22,8 +23,9 @@ subroutine get_obs_sat
!---------------------------------------------------------------------------
!-----------------------------------------------------------------------
- ! !
- ! Load Chlorophyll observations !
+ ! Version 1: A. Teruzzi 2018 !
+ !
+ ! Load Chlorophyll observations !
! !
!-----------------------------------------------------------------------
@@ -33,6 +35,7 @@ subroutine get_obs_sat
use obs_str
use mpi_str
use pnetcdf
+ use da_params
use filenames
implicit none
@@ -41,13 +44,29 @@ subroutine get_obs_sat
INTEGER(i4) :: i
REAL(r8) :: zbo, zbn
REAL(r4), ALLOCATABLE :: chl_mis(:,:),chl_err(:,:)
- INTEGER(i4) :: stat, ncid, idvar, VarId
+ REAL(i4), ALLOCATABLE :: flagblk(:,:)
+ INTEGER(i4) :: stat, ncid, idvar, VarId, ierr
+ INTEGER(i4) :: xid, yid, idF, ii
+ INTEGER(KIND=MPI_OFFSET_KIND) :: global_im, global_jm
+ INTEGER(KIND=MPI_OFFSET_KIND) :: MyCountTwod(2), MyStartTwod(2)
+ CHARACTER(LEN=45) :: flagFile
+ CHARACTER(LEN=15) :: FlagVarName
+
+ global_im = GlobalRow
+ global_jm = GlobalCol
+ do ii=1,2
+ MyCountTwod(ii) = MyCount(ii)
+ MyStartTwod(ii) = MyStart(ii)
+ enddo
+
+
+
sat%no = 0
sat%nc = 0
stat = nf90mpi_open(Var3DCommunicator, trim(MISFIT_FILE), NF90_NOWRITE, MPI_INFO_NULL, ncid)
- if (stat .ne. NF90_NOERR ) call handle_err('nf90mpi_open', stat)
+ if (stat .ne. NF90_NOERR ) call handle_err('nf90mpi_open '//trim(MISFIT_FILE), stat)
if(stat.ne.0)then
sat%no = 0
@@ -76,6 +95,7 @@ subroutine get_obs_sat
ALLOCATE ( chl_mis(grd%im,grd%jm) ) ; chl_mis = huge(chl_mis(1,1))
ALLOCATE ( chl_err(grd%im,grd%jm) ) ; chl_err = huge(chl_err(1,1))
+ ALLOCATE ( flagblk(grd%im,grd%jm) )
ALLOCATE ( sat%flg(sat%no)) ; sat%flg = huge(sat%flg(1))
ALLOCATE ( sat%flc(sat%no)) ; sat%flc = huge(sat%flc(1))
ALLOCATE ( sat%inc(sat%no)) ; sat%inc = huge(sat%inc(1))
@@ -108,6 +128,7 @@ subroutine get_obs_sat
sat%res(k) = chl_mis(i,j)
sat%err(k) = chl_err(i,j)
enddo
+
! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST
! sat%res(:) = 0.
@@ -123,6 +144,7 @@ subroutine get_obs_sat
! ---
! Initialise quality flag, do residual check, compute vertical integration parameters and count good observations
+ flagblk(:,:) = 0 !blacklisting flag
sat%nc = 0
do k=1,sat%no
j = (k-1)/grd%im + 1
@@ -132,6 +154,9 @@ subroutine get_obs_sat
if(abs(sat%res(k)).gt.sat%max_val) then
! residual check
sat%flg(k) = 0
+ if(abs(sat%res(k)).lt.100) then
+ flagblk(i,j) = 1
+ endif
else
! compute vertical integration parameters
zbn = grd%dep(1)*2.0
@@ -151,7 +176,40 @@ subroutine get_obs_sat
if(MyId .eq. 0) then
print*,'Good chl observations: ',sat%nc_global
+ print*,'Saving flag misfit'
endif
+
+ ! Saving flag misfit sat
+ flagFile = 'DA__FREQ_1/flagsat.'//ShortDate//'.nc'
+ ierr = nf90mpi_create(Var3DCommunicator, trim(flagFile), NF90_CLOBBER, MPI_INFO_NULL,ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('flagFile ', ierr)
+
+ ierr = nf90mpi_def_dim(ncid,'x',global_im ,xid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim longitude ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'y' ,global_jm ,yid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim latitude ', ierr)
+
+ FlagVarName='flag_lim_misf'
+
+ ierr = nf90mpi_def_var(ncid, FlagVarName, nf90_int, (/xid,yid/), idF )
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', ierr)
+
+ ierr = nf90mpi_enddef(ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_enddef', ierr)
+
+
+ ierr = nf90mpi_put_var_all(ncid,idF,flagblk,MyStartTwod,MyCountTwod)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all ', ierr)
+
+ ierr = nf90mpi_close(ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('LimCorrfile ', ierr)
+
+
+
+ DEALLOCATE ( flagblk )
+
+
+
sat%flc(:) = sat%flg(:)
end subroutine get_obs_sat
@@ -162,7 +220,7 @@ subroutine int_par_chl
! !
! Get interpolation parameters for a grid !
! !
- ! Version 1: S.Dobricic 2006 !
+ ! Version 1: A. Teruzzi 2018 !
!-----------------------------------------------------------------------
use set_knd
diff --git a/grd_str.f90 b/grd_str.f90
index c977c12..c23b2b1 100644
--- a/grd_str.f90
+++ b/grd_str.f90
@@ -71,6 +71,8 @@ MODULE grd_str
INTEGER(i4) :: imax ! Maximum number of extended points
INTEGER(i4) :: jmax ! Maximum number of extended points
+ INTEGER(i4) :: imax_chl ! Maximum number of extended points until depth pf chl for multivariate
+ INTEGER(i4) :: jmax_chl ! Maximum number of extended points until depth of chl for multivariete
INTEGER(i4), POINTER :: imx(:) ! Max. no. of extended pnts at each level
INTEGER(i4), POINTER :: jmx(:) ! Max. no. of extended pnts at each level
INTEGER(i4), POINTER :: istp(:,:,:) ! Extended points
diff --git a/main.f90 b/main.f90
index 2b5747c..af4dded 100644
--- a/main.f90
+++ b/main.f90
@@ -15,6 +15,7 @@ program ocean_var
call oceanvar
! finalizing the MPI environment
+call clean_da_params
call mpi_stop
end program ocean_var
diff --git a/mpi_str.f90 b/mpi_str.f90
index facc871..7858202 100644
--- a/mpi_str.f90
+++ b/mpi_str.f90
@@ -53,10 +53,142 @@ MODULE mpi_str
! Arrays needed for alltoallv communication
! X dimension
integer, allocatable, dimension(:) :: SendCountX2D, SendCountX3D, SendDisplX2D, SendDisplX3D
+ integer, allocatable, dimension(:) :: SendCountX3D_chl, SendDisplX3D_chl
integer, allocatable, dimension(:) :: RecCountX2D, RecCountX3D, RecDisplX2D, RecDisplX3D
+ integer, allocatable, dimension(:) :: RecCountX3D_chl, RecDisplX3D_chl
! Arrays needed for the ghost cells exchange
REAL(r8), POINTER, DIMENSION(:,:) :: ChlExtended
REAL(r8), POINTER, DIMENSION(:) :: SendTop, RecBottom, SendBottom, RecTop
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: SendTop_2d, RecBottom_2d
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: SendBottom_2d, RecTop_2d
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: ChlExtended_3d, N3nExtended_3d, O2oExtended_3d
+
+
+CONTAINS
+
+ SUBROUTINE EXTEND_2D(INPUT, my_km, OUTPUT_Extended)
+ use set_knd
+ use grd_str
+ use obs_str
+ use drv_str
+ use bio_str
+ IMPLICIT NONE
+ INTEGER(i4) :: my_km
+ REAL(r8), DIMENSION(grd%im , grd%jm , my_km), INTENT(IN) :: INPUT
+ REAL(r8), DIMENSION(grd%im+1, grd%jm , my_km), INTENT(OUT) :: OUTPUT_Extended
+
+
+ INTEGER(i4) :: i, j, k, kk
+ INTEGER :: MyTag
+ INTEGER :: ReqTop, ReqBottom, ierr
+ INTEGER :: StatBottom(MPI_STATUS_SIZE)
+
+
+ ! Filling array to send
+ do k=1,my_km
+ do j=1,grd%jm
+ SendTop_2d(j,k) = INPUT(1,j,k)
+ end do
+ end do
+
+ do k=my_km+1,grd%km
+ SendTop_2d(:,k) = 0
+ end do
+
+ MyTag = 42
+ RecBottom_2d(:,:) = 0
+
+ call MPI_Isend(SendTop_2d, grd%jm*grd%km, MPI_REAL8, ProcTop, MyTag, &
+ Var3DCommunicator, ReqTop, ierr)
+ call MPI_Irecv(RecBottom_2d, grd%jm*grd%km, MPI_REAL8, ProcBottom, MyTag, &
+ Var3DCommunicator, ReqBottom, ierr)
+
+ do k=1,my_km
+ do j=1,grd%jm
+ do i=1,grd%im
+ OUTPUT_Extended(i,j,k) = INPUT(i,j,k)
+ end do
+ end do
+ end do
+
+
+ call MPI_Wait(ReqBottom, StatBottom, ierr)
+ do k=1,my_km
+ do j=1,grd%jm
+ OUTPUT_Extended(grd%im+1,j,k) = RecBottom_2d(j,k)
+ end do
+ end do
+
+
+
+ END SUBROUTINE EXTEND_2D
+
+
+ SUBROUTINE ADD_PREVCORE_CONTRIB(INPUT_Extended,my_km,OUTPUT,INIT_2d)
+ ! ADDS PREVIOUS CORE CONTRIBUTION to the sum of an _ad variable
+ ! used in obs_arg_ad
+ !
+ ! THEORY
+
+ ! OUTPUT = INIT + contr(i) + contr(i+1)
+ ! but without ghost cell we have
+ ! OUTPUT(1) = INIT + contr(1)
+ ! OUTPUT(im+1) = INIT + contr(im+1)
+ ! i.e one single contribution
+
+ ! im+1 is position 1 of following core,
+ ! then we can receive contr(im+1) in RecTop_2d
+
+ ! In order to have OUTPUT = INIT + contr(i) + contr(i+1) we'll do
+ ! OUTPUT(1) = INIT + contrib(1) + INIT + contr(i+1) - INIT
+ ! OUTPUT(1) + Rec_top - INIT
+
+ use set_knd
+ use grd_str
+ use obs_str
+ use drv_str
+ use bio_str
+ IMPLICIT NONE
+ INTEGER(i4) :: my_km
+ REAL(r8), DIMENSION(grd%im+1, grd%jm , my_km), INTENT(IN) :: INPUT_Extended
+ REAL(r8), DIMENSION(grd%im , grd%jm , my_km), INTENT(OUT) :: OUTPUT
+ REAL(r8), DIMENSION( grd%jm , my_km), INTENT(IN) :: INIT_2d
+
+ INTEGER :: ReqBottom, ReqTop, ierr
+ INTEGER(i4) :: i, j, k, kk
+ INTEGER :: MyTag
+ INTEGER :: StatTop(MPI_STATUS_SIZE)
+
+ do k=1,my_km
+ do j=1,grd%jm
+ SendBottom_2d(j,k) = INPUT_Extended(grd%im+1,j,k)
+ end do
+ end do
+
+
+ do k=my_km+1,grd%km
+ SendBottom_2d(:,k) = 0
+ end do
+
+ MyTag = 42
+ RecTop_2d(:,:) = 0
+
+ call MPI_Isend(SendBottom_2d, grd%jm*grd%km, MPI_REAL8, ProcBottom, MyTag, &
+ Var3DCommunicator, ReqBottom, ierr)
+ call MPI_Irecv(RecTop_2d, grd%jm*grd%km, MPI_REAL8, ProcTop, MyTag, &
+ Var3DCommunicator, ReqTop, ierr)
+
+ OUTPUT = INPUT_Extended(1:grd%im,:,:)
+
+ call MPI_Wait(ReqTop, StatTop, ierr)
+
+ do k=1,my_km
+ do j=1,grd%jm
+ OUTPUT(1,j,k) = OUTPUT(1,j,k) + RecTop_2d(j,k) - INIT_2d(j,k)
+ end do
+ end do
+
+ END SUBROUTINE ADD_PREVCORE_CONTRIB
END MODULE mpi_str
diff --git a/mpi_utils.f90 b/mpi_utils.f90
index 17ff534..443c9f0 100644
--- a/mpi_utils.f90
+++ b/mpi_utils.f90
@@ -1,6 +1,12 @@
subroutine var3d_mpi_init()
+#include
+#if PETSC_VERSION_GE(3,8,0)
+#include "petsc/finclude/petscvec.h"
+#else
#include "petsc/finclude/petscvecdef.h"
+#endif
+
use mpi_str
use drv_str
use petscvec
@@ -93,9 +99,13 @@ subroutine my_3dvar_node()
call MPI_TYPE_COMMIT(MyPair, ierr)
ALLOCATE(SendCountX2D(NumProcI), SendCountX3D(NumProcI))
+ ALLOCATE(SendCountX3D_chl(NumProcI))
ALLOCATE(SendDisplX2D(NumProcI), SendDisplX3D(NumProcI))
+ ALLOCATE(SendDisplX3D_chl(NumProcI))
ALLOCATE(RecCountX2D(NumProcI), RecCountX3D(NumProcI))
+ ALLOCATE(RecCountX3D_chl(NumProcI))
ALLOCATE(RecDisplX2D(NumProcI), RecDisplX3D(NumProcI))
+ ALLOCATE(RecDisplX3D_chl(NumProcI))
! print for debug
! write(*,*) "MyId", MyId, "PosI", MyPosI, "PosJ", MyPosJ, "Left", ProcLeft, "Right", ProcRight, "Top", ProcTop, "Bottom", ProcBottom
@@ -136,7 +146,13 @@ end subroutine mpi_sync
subroutine mpi_stop
+#include
+#if PETSC_VERSION_GE(3,8,0)
+#include "petsc/finclude/petscvec.h"
+#else
#include "petsc/finclude/petscvecdef.h"
+#endif
+
use mpi_str
use petscvec
diff --git a/namelists/satfloat.20150101.nml b/namelists/satfloat.20150101.nml
new file mode 100644
index 0000000..5fbd536
--- /dev/null
+++ b/namelists/satfloat.20150101.nml
@@ -0,0 +1,59 @@
+!------------------------------------------------------------
+! OceanVar namelist specification for multivariate and multiplatform
+!------------------------------------------------------------
+!------------------------------------------------------------
+!
+! Namelist biolst
+! ---
+!
+! Biological assimilation set-up
+!
+! chl_assim - Chlorophyll assimilation
+! chl_upnut - Nutrient update based on chl assimilation chl_upnut
+! nut - Nutrient assimilation
+! multiv - Multivariate assimilation
+! nphyto - Number of phytoplankton species
+! chl_dep - Minimum depth for chlorophyll
+! ncmp - Number of phytoplankton components
+! ApplyConditions - Apply conditions flag
+! N3n - N3n assimilation
+! updateN1p - N1p update based on N3n assimilation updateN1p
+! O2o - O2o assimilation
+!
+! ---
+&biolst
+ chl_assim = 1
+ chl_upnut = 0
+ nut = 1
+ nphyto = 4
+ chl_dep = 0.0
+ ncmp = 5
+ApplyConditions = .true.
+ N3n = 1
+ updateN1p = 0
+ O2o = 1
+/
+!------------------------------------------------------------
+!
+! Namelist parameters
+! ---
+!
+! Parameters namelist
+!
+! sat - 1-assimilate satellite data
+! 0-no satellite assimilation
+! argo - 1-assimilate argo data
+! - 0-no argo assimilation
+! uniformL - 1-non uniform radius
+! - 0-uniform radius (rcf%L)
+! anisL - 1-anisotropy
+! 0-isotropy
+!
+! ---
+¶ms
+ sat_obs = 1
+ argo = 1
+ uniformL = 0
+ anisL = 0
+/
+!------------------------------------------------------------
diff --git a/namelists/var_3d_nml b/namelists/var_3d_nml
index 40f756b..a2cbc7d 100644
--- a/namelists/var_3d_nml
+++ b/namelists/var_3d_nml
@@ -8,14 +8,15 @@
!
! BFGS minimizers set-up
!
-! ctl_m - Number of copies saved in the minimizer
-! ctl_tol - Stopping criteria (absolute)
-! ctl_per - Stopping criteria (relative)
+! ctl_tol - Minimum gradient of J
+! ctl_per - Percentage of initial gradient
+! verbose - Add verbose on standard output
!
! ---
&ctllst
ctl_tol = 1.e-11
ctl_per = 1.e-3
+ verbose = 1
/
!------------------------------------------------------------
!
@@ -24,8 +25,13 @@
!
! Covariance constants
!
-! neof - Number of vertical EOFs
-! nreg - Number of regions
+! neof_chl - Number of vertical EOFs for chl
+! neof_n3n - Number of vertical EOFs for n3n
+! neof_o2o - Number of vertical EOFs for o2o
+! neof_multi - Number of multivariate vertical EOFs
+! kmchl - chl Nlevels in multi EOFs
+! kmnit - nit Nlevels in multi EOFs
+! nreg - Number of regions
! read_eof - Logical to read EOFs from files.
! See subroutine def_cov.f90
! rcf_ntr - Number of iterations of the recursive filter
@@ -37,64 +43,13 @@
neof_chl = 4
neof_n3n = 4
neof_o2o = 4
- nreg = 63045
+ neof_multi = 26
+ kmchl = 26
+ kmnit = 40
+ nreg = 15
read_eof = .true.
rcf_ntr = 4
rcf_L = 5000.
rcf_efc = 5.0
/
!------------------------------------------------------------
-!
-! Namelist biolst
-! ---
-!
-! Biological assimilation set-up
-!
-! chl_assim - Chlorophyll assimilation
-! nut - Nutrient assimilation
-! nphyto - Number of phytoplankton species
-! chl_dep - Minimum depth for chlorophyll
-! ncmp - Number of phytoplankton components
-! ApplyConditions - Apply conditions flag
-! N3n - N3n assimilation
-! updateN1p - N1p update based on N3n assimilation updateN1p
-! O2o - O2o assimilation
-!
-! ---
-&biolst
- chl_assim = 1
- nut = 1
- nphyto = 4
- chl_dep = 0.0
- ncmp = 5
-ApplyConditions = .true.
- N3n = 1
- updateN1p = 0
- O2o = 1
-/
-!------------------------------------------------------------
-!
-! Namelist parameters
-! ---
-!
-! Parameters namelist
-!
-! sat - 1-assimilate satellite data
-! 0-no satellite assimilation
-! argo - 1-assimilate argo data
-! - 0-no argo assimilation
-! uniformL - 1-non uniform radius
-! - 0-uniform radius (rcf%L)
-! anisL - 1-anisotropy
-! 0-isotropy
-! verbose - 1-set verbose output
-!
-! ---
-¶ms
- sat_obs = 1
- argo = 1
- uniformL = 0
- anisL = 0
- verbose = 1
-/
-!------------------------------------------------------------
diff --git a/obs_arg.f90 b/obs_arg.f90
index 49e3f09..39d95f8 100644
--- a/obs_arg.f90
+++ b/obs_arg.f90
@@ -5,7 +5,7 @@ subroutine obs_arg
! !
! Copyright 2006 Srdjan Dobricic, CMCC, Bologna !
! !
- ! This file is part of OceanVar. !
+ ! This file is part of OceanVar. !
! !
! OceanVar is free software: you can redistribute it and/or modify. !
! it under the terms of the GNU General Public License as published by !
@@ -18,7 +18,7 @@ subroutine obs_arg
! GNU General Public License for more details. !
! !
! You should have received a copy of the GNU General Public License !
- ! along with OceanVar. If not, see . !
+ ! along with OceanVar. If not, see . !
! !
!---------------------------------------------------------------------------
@@ -32,137 +32,86 @@ subroutine obs_arg
use set_knd
use grd_str
+ use eof_str
use obs_str
use mpi_str
use drv_str
use bio_str
implicit none
+
+ INTEGER(i4) :: i, j, k, kk, condc, condn, my_km
- INTEGER(i4) :: i, j, k, kk
- REAL(r8) :: FirstData, SecData, ThirdData, LastData, Test
- INTEGER(kind=MPI_ADDRESS_KIND) :: TargetOffset
- INTEGER :: ierr, NData
- real(r8), pointer, dimension(:,:,:) :: GetData
-
+ my_km = grd%km
+ if(drv%multiv.eq.1) &
+ my_km = ros%kmchl
+
+ condc = 0
+ condn = 0
+ if ((drv%chl_assim.eq.1 ) .or. (drv%multiv.eq.1)) then
+ condc = 1
+ call EXTEND_2D( grd%chl, my_km, ChlExtended_3d )
+ endif
+ if ((drv%nut.eq.1 .and. bio%N3n.eq.1 ) .or. (drv%multiv.eq.1)) then
+ condn = 1
+ call EXTEND_2D( grd%n3n, grd%km, N3nExtended_3d )
+ endif
+ if (bio%O2o.eq.1 ) &
+ call EXTEND_2D( grd%O2o, grd%km, O2oExtended_3d )
+
+
+
do kk = 1,arg%no
-
- if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.0 .and. drv%chl_assim.eq.1)then
-
- i=arg%ib(kk)
- j=arg%jb(kk)
- k=arg%kb(kk)
-
- if(i .lt. grd%im) then
- arg%inc(kk) = arg%pq1(kk) * grd%chl(i ,j ,k) + &
- arg%pq2(kk) * grd%chl(i+1,j ,k ) + &
- arg%pq3(kk) * grd%chl(i ,j+1,k ) + &
- arg%pq4(kk) * grd%chl(i+1,j+1,k ) + &
- arg%pq5(kk) * grd%chl(i ,j ,k+1) + &
- arg%pq6(kk) * grd%chl(i+1,j ,k+1) + &
- arg%pq7(kk) * grd%chl(i ,j+1,k+1) + &
- arg%pq8(kk) * grd%chl(i+1,j+1,k+1)
-
- else
- ALLOCATE(GetData(NextLocalRow,grd%jm,2))
-
- NData = NextLocalRow*grd%jm*2
- call MPI_Win_lock (MPI_LOCK_EXCLUSIVE, ProcBottom, 0, MpiWinChl, ierr )
- TargetOffset = (k-1)*grd%jm*NextLocalRow
- call MPI_Get (GetData, NData, MPI_REAL8, ProcBottom, TargetOffset, NData, MPI_REAL8, MpiWinChl, ierr)
- call MPI_Win_unlock(ProcBottom, MpiWinChl, ierr)
-
- arg%inc(kk) = arg%pq1(kk) * grd%chl(i ,j ,k) + &
- arg%pq2(kk) * GetData(1 ,j ,1 ) + &
- arg%pq3(kk) * grd%chl(i ,j+1,k ) + &
- arg%pq4(kk) * GetData(1 ,j+1,1 ) + &
- arg%pq5(kk) * grd%chl(i ,j ,k+1) + &
- arg%pq6(kk) * GetData(1 ,j ,2 ) + &
- arg%pq7(kk) * grd%chl(i ,j+1,k+1) + &
- arg%pq8(kk) * GetData(1 ,j+1,2 )
-
-
- DEALLOCATE(GetData)
- endif
-
- else if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.1 .and. drv%nut.eq.1 .and. bio%n3n.eq.1) then
-
- i=arg%ib(kk)
- j=arg%jb(kk)
- k=arg%kb(kk)
-
- if(i .lt. grd%im) then
- arg%inc(kk) = arg%pq1(kk) * grd%n3n(i ,j ,k) + &
- arg%pq2(kk) * grd%n3n(i+1,j ,k ) + &
- arg%pq3(kk) * grd%n3n(i ,j+1,k ) + &
- arg%pq4(kk) * grd%n3n(i+1,j+1,k ) + &
- arg%pq5(kk) * grd%n3n(i ,j ,k+1) + &
- arg%pq6(kk) * grd%n3n(i+1,j ,k+1) + &
- arg%pq7(kk) * grd%n3n(i ,j+1,k+1) + &
- arg%pq8(kk) * grd%n3n(i+1,j+1,k+1)
-
- else
- ALLOCATE(GetData(NextLocalRow,grd%jm,2))
-
- NData = NextLocalRow*grd%jm*2
- call MPI_Win_lock (MPI_LOCK_EXCLUSIVE, ProcBottom, 0, MpiWinN3n, ierr )
- TargetOffset = (k-1)*grd%jm*NextLocalRow
- call MPI_Get (GetData, NData, MPI_REAL8, ProcBottom, TargetOffset, NData, MPI_REAL8, MpiWinN3n, ierr)
- call MPI_Win_unlock(ProcBottom, MpiWinN3n, ierr)
-
- arg%inc(kk) = arg%pq1(kk) * grd%n3n(i ,j ,k) + &
- arg%pq2(kk) * GetData(1 ,j ,1 ) + &
- arg%pq3(kk) * grd%n3n(i ,j+1,k ) + &
- arg%pq4(kk) * GetData(1 ,j+1,1 ) + &
- arg%pq5(kk) * grd%n3n(i ,j ,k+1) + &
- arg%pq6(kk) * GetData(1 ,j ,2 ) + &
- arg%pq7(kk) * grd%n3n(i ,j+1,k+1) + &
- arg%pq8(kk) * GetData(1 ,j+1,2 )
-
-
- DEALLOCATE(GetData)
- endif
-
- else if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.2 .and. drv%nut.eq.1 .and. bio%o2o.eq.1) then
-
+
i=arg%ib(kk)
j=arg%jb(kk)
k=arg%kb(kk)
-
- if(i .lt. grd%im) then
- arg%inc(kk) = arg%pq1(kk) * grd%o2o(i ,j ,k) + &
- arg%pq2(kk) * grd%o2o(i+1,j ,k ) + &
- arg%pq3(kk) * grd%o2o(i ,j+1,k ) + &
- arg%pq4(kk) * grd%o2o(i+1,j+1,k ) + &
- arg%pq5(kk) * grd%o2o(i ,j ,k+1) + &
- arg%pq6(kk) * grd%o2o(i+1,j ,k+1) + &
- arg%pq7(kk) * grd%o2o(i ,j+1,k+1) + &
- arg%pq8(kk) * grd%o2o(i+1,j+1,k+1)
-
- else
- ALLOCATE(GetData(NextLocalRow,grd%jm,2))
-
- NData = NextLocalRow*grd%jm*2
- call MPI_Win_lock (MPI_LOCK_EXCLUSIVE, ProcBottom, 0, MpiWinO2o, ierr )
- TargetOffset = (k-1)*grd%jm*NextLocalRow
- call MPI_Get (GetData, NData, MPI_REAL8, ProcBottom, TargetOffset, NData, MPI_REAL8, MpiWinO2o, ierr)
- call MPI_Win_unlock(ProcBottom, MpiWinO2o, ierr)
-
- arg%inc(kk) = arg%pq1(kk) * grd%o2o(i ,j ,k) + &
- arg%pq2(kk) * GetData(1 ,j ,1 ) + &
- arg%pq3(kk) * grd%o2o(i ,j+1,k ) + &
- arg%pq4(kk) * GetData(1 ,j+1,1 ) + &
- arg%pq5(kk) * grd%o2o(i ,j ,k+1) + &
- arg%pq6(kk) * GetData(1 ,j ,2 ) + &
- arg%pq7(kk) * grd%o2o(i ,j+1,k+1) + &
- arg%pq8(kk) * GetData(1 ,j+1,2 )
-
-
- DEALLOCATE(GetData)
- endif
-
+
+ if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.0 .and. condc.eq.1) then
+
+ arg%inc(kk) = &
+ arg%pq1(kk) * ChlExtended_3d(i ,j ,k) + &
+ arg%pq2(kk) * ChlExtended_3d(i+1,j ,k ) + &
+ arg%pq3(kk) * ChlExtended_3d(i ,j+1,k ) + &
+ arg%pq4(kk) * ChlExtended_3d(i+1,j+1,k ) + &
+ arg%pq5(kk) * ChlExtended_3d(i ,j ,k+1) + &
+ arg%pq6(kk) * ChlExtended_3d(i+1,j ,k+1) + &
+ arg%pq7(kk) * ChlExtended_3d(i ,j+1,k+1) + &
+ arg%pq8(kk) * ChlExtended_3d(i+1,j+1,k+1)
+ ! if(drv%multiv .eq. 1) &
+ ! arg%inc(kk) = arg%inc(kk) * arg%std(kk)
+ endif
+ if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.1 .and. condn.eq.1) then
+
+
+ arg%inc(kk) = &
+ arg%pq1(kk) * N3nExtended_3d(i ,j ,k) + &
+ arg%pq2(kk) * N3nExtended_3d(i+1,j ,k ) + &
+ arg%pq3(kk) * N3nExtended_3d(i ,j+1,k ) + &
+ arg%pq4(kk) * N3nExtended_3d(i+1,j+1,k ) + &
+ arg%pq5(kk) * N3nExtended_3d(i ,j ,k+1) + &
+ arg%pq6(kk) * N3nExtended_3d(i+1,j ,k+1) + &
+ arg%pq7(kk) * N3nExtended_3d(i ,j+1,k+1) + &
+ arg%pq8(kk) * N3nExtended_3d(i+1,j+1,k+1)
+ ! if(drv%multiv .eq. 1) &
+ ! arg%inc(kk) = arg%inc(kk) * arg%std(kk)
endif
-
+ if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.2 .and. drv%nut.eq.1 .and. bio%o2o.eq.1) then
+
+ arg%inc(kk) = &
+ arg%pq1(kk) * O2oExtended_3d(i ,j ,k) + &
+ arg%pq2(kk) * O2oExtended_3d(i+1,j ,k ) + &
+ arg%pq3(kk) * O2oExtended_3d(i ,j+1,k ) + &
+ arg%pq4(kk) * O2oExtended_3d(i+1,j+1,k ) + &
+ arg%pq5(kk) * O2oExtended_3d(i ,j ,k+1) + &
+ arg%pq6(kk) * O2oExtended_3d(i+1,j ,k+1) + &
+ arg%pq7(kk) * O2oExtended_3d(i ,j+1,k+1) + &
+ arg%pq8(kk) * O2oExtended_3d(i+1,j+1,k+1)
+
+ endif
+
+
enddo
+
end subroutine obs_arg
diff --git a/obs_arg_ad.f90 b/obs_arg_ad.f90
index f1ca211..db4b78d 100644
--- a/obs_arg_ad.f90
+++ b/obs_arg_ad.f90
@@ -5,7 +5,7 @@ subroutine obs_arg_ad
! !
! Copyright 2006 Srdjan Dobricic, CMCC, Bologna !
! !
- ! This file is part of OceanVar. !
+ ! This file is part of OceanVar. !
! !
! OceanVar is free software: you can redistribute it and/or modify. !
! it under the terms of the GNU General Public License as published by !
@@ -18,7 +18,7 @@ subroutine obs_arg_ad
! GNU General Public License for more details. !
! !
! You should have received a copy of the GNU General Public License !
- ! along with OceanVar. If not, see . !
+ ! along with OceanVar. If not, see . !
! !
!---------------------------------------------------------------------------
@@ -32,6 +32,7 @@ subroutine obs_arg_ad
use set_knd
use grd_str
+ use eof_str
use obs_str
use mpi_str
use filenames
@@ -40,148 +41,107 @@ subroutine obs_arg_ad
implicit none
- INTEGER(i4) :: i, j, k, kk
- REAL(r8) :: ToSum
- INTEGER(kind=MPI_ADDRESS_KIND) :: TargetOffset
- INTEGER :: ierr, NData
- real(r8), pointer, dimension(:,:,:) :: MatrixToSum
+ INTEGER(i4) :: i, j, k, kk, condc, condn, my_km
+ REAL(r8), DIMENSION(grd%jm,grd%km) :: slicevar
+ REAL(8) :: obsg
+
+ my_km = grd%km
+ if(drv%multiv.eq.1) &
+ my_km = ros%kmchl
+ condc = 0
+ condn = 0
+ if ((drv%chl_assim.eq.1 ) .or. (drv%multiv.eq.1)) then
+ condc = 1
+ call EXTEND_2D( grd%chl_ad, my_km, ChlExtended_3d )
+ endif
+ if ((drv%nut.eq.1 .and. bio%N3n.eq.1 ) .or. (drv%multiv.eq.1)) then
+ call EXTEND_2D( grd%n3n_ad, grd%km, N3nExtended_3d )
+ condn = 1
+ endif
+ if (drv%nut.eq.1 .and. bio%O2o.eq.1 ) &
+ call EXTEND_2D( grd%O2o_ad, grd%km, O2oExtended_3d )
+
+
do kk = 1,arg%no
-
- if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.0 .and. drv%chl_assim.eq.1)then
-
- obs%k = obs%k + 1
-
- i=arg%ib(kk)
- j=arg%jb(kk)
- k=arg%kb(kk)
-
- if(i .lt. grd%im) then
-
- grd%chl_ad(i ,j ,k ) = grd%chl_ad(i ,j ,k ) + arg%pq1(kk) * obs%gra(obs%k)
- grd%chl_ad(i+1,j ,k ) = grd%chl_ad(i+1,j ,k ) + arg%pq2(kk) * obs%gra(obs%k)
- grd%chl_ad(i ,j+1,k ) = grd%chl_ad(i ,j+1,k ) + arg%pq3(kk) * obs%gra(obs%k)
- grd%chl_ad(i+1,j+1,k ) = grd%chl_ad(i+1,j+1,k ) + arg%pq4(kk) * obs%gra(obs%k)
- grd%chl_ad(i ,j ,k+1) = grd%chl_ad(i ,j ,k+1) + arg%pq5(kk) * obs%gra(obs%k)
- grd%chl_ad(i+1,j ,k+1) = grd%chl_ad(i+1,j ,k+1) + arg%pq6(kk) * obs%gra(obs%k)
- grd%chl_ad(i ,j+1,k+1) = grd%chl_ad(i ,j+1,k+1) + arg%pq7(kk) * obs%gra(obs%k)
- grd%chl_ad(i+1,j+1,k+1) = grd%chl_ad(i+1,j+1,k+1) + arg%pq8(kk) * obs%gra(obs%k)
-
- else
-
- ALLOCATE(MatrixToSum(NextLocalRow,grd%jm,2))
- MatrixToSum(:,:,:) = dble(0)
-
- grd%chl_ad(i ,j ,k ) = grd%chl_ad(i ,j ,k ) + arg%pq1(kk) * obs%gra(obs%k)
- grd%chl_ad(i ,j+1,k ) = grd%chl_ad(i ,j+1,k ) + arg%pq3(kk) * obs%gra(obs%k)
- grd%chl_ad(i ,j ,k+1) = grd%chl_ad(i ,j ,k+1) + arg%pq5(kk) * obs%gra(obs%k)
- grd%chl_ad(i ,j+1,k+1) = grd%chl_ad(i ,j+1,k+1) + arg%pq7(kk) * obs%gra(obs%k)
-
- MatrixToSum(1,j ,1) = arg%pq2(kk) * obs%gra(obs%k)
- MatrixToSum(1,j+1,1) = arg%pq4(kk) * obs%gra(obs%k)
- MatrixToSum(1,j ,2) = arg%pq6(kk) * obs%gra(obs%k)
- MatrixToSum(1,j+1,2) = arg%pq8(kk) * obs%gra(obs%k)
-
- call MPI_Win_lock (MPI_LOCK_EXCLUSIVE, ProcBottom, 0, MpiWinChlAd, ierr )
- NData = NextLocalRow*grd%jm*2
- TargetOffset = (k-1)*grd%jm*NextLocalRow
- call MPI_Accumulate (MatrixToSum, NData, MPI_REAL8, ProcBottom, TargetOffset, NData, MPI_REAL8, MPI_SUM, MpiWinChlAd, ierr)
-
- call MPI_Win_unlock(ProcBottom, MpiWinChlAd, ierr)
- DEALLOCATE(MatrixToSum)
-
- endif
-
- else if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.1 .and. drv%nut.eq.1 .and. bio%n3n.eq.1) then
-
- obs%k = obs%k + 1
-
- i=arg%ib(kk)
- j=arg%jb(kk)
- k=arg%kb(kk)
-
- if(i .lt. grd%im) then
-
- grd%n3n_ad(i ,j ,k ) = grd%n3n_ad(i ,j ,k ) + arg%pq1(kk) * obs%gra(obs%k)
- grd%n3n_ad(i+1,j ,k ) = grd%n3n_ad(i+1,j ,k ) + arg%pq2(kk) * obs%gra(obs%k)
- grd%n3n_ad(i ,j+1,k ) = grd%n3n_ad(i ,j+1,k ) + arg%pq3(kk) * obs%gra(obs%k)
- grd%n3n_ad(i+1,j+1,k ) = grd%n3n_ad(i+1,j+1,k ) + arg%pq4(kk) * obs%gra(obs%k)
- grd%n3n_ad(i ,j ,k+1) = grd%n3n_ad(i ,j ,k+1) + arg%pq5(kk) * obs%gra(obs%k)
- grd%n3n_ad(i+1,j ,k+1) = grd%n3n_ad(i+1,j ,k+1) + arg%pq6(kk) * obs%gra(obs%k)
- grd%n3n_ad(i ,j+1,k+1) = grd%n3n_ad(i ,j+1,k+1) + arg%pq7(kk) * obs%gra(obs%k)
- grd%n3n_ad(i+1,j+1,k+1) = grd%n3n_ad(i+1,j+1,k+1) + arg%pq8(kk) * obs%gra(obs%k)
-
- else
-
- ALLOCATE(MatrixToSum(NextLocalRow,grd%jm,2))
- MatrixToSum(:,:,:) = dble(0)
-
- grd%n3n_ad(i ,j ,k ) = grd%n3n_ad(i ,j ,k ) + arg%pq1(kk) * obs%gra(obs%k)
- grd%n3n_ad(i ,j+1,k ) = grd%n3n_ad(i ,j+1,k ) + arg%pq3(kk) * obs%gra(obs%k)
- grd%n3n_ad(i ,j ,k+1) = grd%n3n_ad(i ,j ,k+1) + arg%pq5(kk) * obs%gra(obs%k)
- grd%n3n_ad(i ,j+1,k+1) = grd%n3n_ad(i ,j+1,k+1) + arg%pq7(kk) * obs%gra(obs%k)
-
- MatrixToSum(1,j ,1) = arg%pq2(kk) * obs%gra(obs%k)
- MatrixToSum(1,j+1,1) = arg%pq4(kk) * obs%gra(obs%k)
- MatrixToSum(1,j ,2) = arg%pq6(kk) * obs%gra(obs%k)
- MatrixToSum(1,j+1,2) = arg%pq8(kk) * obs%gra(obs%k)
-
- call MPI_Win_lock (MPI_LOCK_EXCLUSIVE, ProcBottom, 0, MpiWinN3nAd, ierr )
- NData = NextLocalRow*grd%jm*2
- TargetOffset = (k-1)*grd%jm*NextLocalRow
- call MPI_Accumulate (MatrixToSum, NData, MPI_REAL8, ProcBottom, TargetOffset, NData, MPI_REAL8, MPI_SUM, MpiWinN3nAd, ierr)
-
- call MPI_Win_unlock(ProcBottom, MpiWinN3nAd, ierr)
- DEALLOCATE(MatrixToSum)
-
- endif
-
- else if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.2 .and. drv%nut.eq.1 .and. bio%o2o.eq.1) then
-
- obs%k = obs%k + 1
-
- i=arg%ib(kk)
- j=arg%jb(kk)
- k=arg%kb(kk)
-
- if(i .lt. grd%im) then
-
- grd%o2o_ad(i ,j ,k ) = grd%o2o_ad(i ,j ,k ) + arg%pq1(kk) * obs%gra(obs%k)
- grd%o2o_ad(i+1,j ,k ) = grd%o2o_ad(i+1,j ,k ) + arg%pq2(kk) * obs%gra(obs%k)
- grd%o2o_ad(i ,j+1,k ) = grd%o2o_ad(i ,j+1,k ) + arg%pq3(kk) * obs%gra(obs%k)
- grd%o2o_ad(i+1,j+1,k ) = grd%o2o_ad(i+1,j+1,k ) + arg%pq4(kk) * obs%gra(obs%k)
- grd%o2o_ad(i ,j ,k+1) = grd%o2o_ad(i ,j ,k+1) + arg%pq5(kk) * obs%gra(obs%k)
- grd%o2o_ad(i+1,j ,k+1) = grd%o2o_ad(i+1,j ,k+1) + arg%pq6(kk) * obs%gra(obs%k)
- grd%o2o_ad(i ,j+1,k+1) = grd%o2o_ad(i ,j+1,k+1) + arg%pq7(kk) * obs%gra(obs%k)
- grd%o2o_ad(i+1,j+1,k+1) = grd%o2o_ad(i+1,j+1,k+1) + arg%pq8(kk) * obs%gra(obs%k)
-
- else
-
- ALLOCATE(MatrixToSum(NextLocalRow,grd%jm,2))
- MatrixToSum(:,:,:) = dble(0)
-
- grd%o2o_ad(i ,j ,k ) = grd%o2o_ad(i ,j ,k ) + arg%pq1(kk) * obs%gra(obs%k)
- grd%o2o_ad(i ,j+1,k ) = grd%o2o_ad(i ,j+1,k ) + arg%pq3(kk) * obs%gra(obs%k)
- grd%o2o_ad(i ,j ,k+1) = grd%o2o_ad(i ,j ,k+1) + arg%pq5(kk) * obs%gra(obs%k)
- grd%o2o_ad(i ,j+1,k+1) = grd%o2o_ad(i ,j+1,k+1) + arg%pq7(kk) * obs%gra(obs%k)
-
- MatrixToSum(1,j ,1) = arg%pq2(kk) * obs%gra(obs%k)
- MatrixToSum(1,j+1,1) = arg%pq4(kk) * obs%gra(obs%k)
- MatrixToSum(1,j ,2) = arg%pq6(kk) * obs%gra(obs%k)
- MatrixToSum(1,j+1,2) = arg%pq8(kk) * obs%gra(obs%k)
-
- call MPI_Win_lock (MPI_LOCK_EXCLUSIVE, ProcBottom, 0, MpiWinO2oAd, ierr )
- NData = NextLocalRow*grd%jm*2
- TargetOffset = (k-1)*grd%jm*NextLocalRow
- call MPI_Accumulate (MatrixToSum, NData, MPI_REAL8, ProcBottom, TargetOffset, NData, MPI_REAL8, MPI_SUM, MpiWinO2oAd, ierr)
-
- call MPI_Win_unlock(ProcBottom, MpiWinO2oAd, ierr)
- DEALLOCATE(MatrixToSum)
-
- endif
-
+
+ i=arg%ib(kk)
+ j=arg%jb(kk)
+ k=arg%kb(kk)
+
+ if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.0 .and. condc.eq.1)then
+
+ obs%k = obs%k + 1
+ ! if(drv%multiv .eq. 1) then
+ ! obsg = obs%gra(obs%k)*arg%std(kk)
+ ! else
+ obsg = obs%gra(obs%k)
+ ! end if
+
+ ChlExtended_3d(i ,j ,k ) = ChlExtended_3d(i ,j ,k ) + arg%pq1(kk) * obsg
+ ChlExtended_3d(i+1,j ,k ) = ChlExtended_3d(i+1,j ,k ) + arg%pq2(kk) * obsg
+ ChlExtended_3d(i ,j+1,k ) = ChlExtended_3d(i ,j+1,k ) + arg%pq3(kk) * obsg
+ ChlExtended_3d(i+1,j+1,k ) = ChlExtended_3d(i+1,j+1,k ) + arg%pq4(kk) * obsg
+ ChlExtended_3d(i ,j ,k+1) = ChlExtended_3d(i ,j ,k+1) + arg%pq5(kk) * obsg
+ ChlExtended_3d(i+1,j ,k+1) = ChlExtended_3d(i+1,j ,k+1) + arg%pq6(kk) * obsg
+ ChlExtended_3d(i ,j+1,k+1) = ChlExtended_3d(i ,j+1,k+1) + arg%pq7(kk) * obsg
+ ChlExtended_3d(i+1,j+1,k+1) = ChlExtended_3d(i+1,j+1,k+1) + arg%pq8(kk) * obsg
+ endif
+ if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.1 .and. condn.eq.1) then
+
+ obs%k = obs%k + 1
+ ! if(drv%multiv .eq. 1) then
+ ! obsg = obs%gra(obs%k)*arg%std(kk)
+ ! else
+ obsg = obs%gra(obs%k)
+ ! end if
+
+ N3nExtended_3d(i ,j ,k ) = N3nExtended_3d(i ,j ,k ) + arg%pq1(kk) * obsg
+ N3nExtended_3d(i+1,j ,k ) = N3nExtended_3d(i+1,j ,k ) + arg%pq2(kk) * obsg
+ N3nExtended_3d(i ,j+1,k ) = N3nExtended_3d(i ,j+1,k ) + arg%pq3(kk) * obsg
+ N3nExtended_3d(i+1,j+1,k ) = N3nExtended_3d(i+1,j+1,k ) + arg%pq4(kk) * obsg
+ N3nExtended_3d(i ,j ,k+1) = N3nExtended_3d(i ,j ,k+1) + arg%pq5(kk) * obsg
+ N3nExtended_3d(i+1,j ,k+1) = N3nExtended_3d(i+1,j ,k+1) + arg%pq6(kk) * obsg
+ N3nExtended_3d(i ,j+1,k+1) = N3nExtended_3d(i ,j+1,k+1) + arg%pq7(kk) * obsg
+ N3nExtended_3d(i+1,j+1,k+1) = N3nExtended_3d(i+1,j+1,k+1) + arg%pq8(kk) * obsg
+
+
endif
-
+ if(arg%flc(kk).eq.1 .and. arg%par(kk).eq.2 .and. drv%nut.eq.1 .and. bio%o2o.eq.1) then
+
+ obs%k = obs%k + 1
+
+ O2oExtended_3d(i ,j ,k ) = O2oExtended_3d(i ,j ,k ) + arg%pq1(kk) * obs%gra(obs%k)
+ O2oExtended_3d(i+1,j ,k ) = O2oExtended_3d(i+1,j ,k ) + arg%pq2(kk) * obs%gra(obs%k)
+ O2oExtended_3d(i ,j+1,k ) = O2oExtended_3d(i ,j+1,k ) + arg%pq3(kk) * obs%gra(obs%k)
+ O2oExtended_3d(i+1,j+1,k ) = O2oExtended_3d(i+1,j+1,k ) + arg%pq4(kk) * obs%gra(obs%k)
+ O2oExtended_3d(i ,j ,k+1) = O2oExtended_3d(i ,j ,k+1) + arg%pq5(kk) * obs%gra(obs%k)
+ O2oExtended_3d(i+1,j ,k+1) = O2oExtended_3d(i+1,j ,k+1) + arg%pq6(kk) * obs%gra(obs%k)
+ O2oExtended_3d(i ,j+1,k+1) = O2oExtended_3d(i ,j+1,k+1) + arg%pq7(kk) * obs%gra(obs%k)
+ O2oExtended_3d(i+1,j+1,k+1) = O2oExtended_3d(i+1,j+1,k+1) + arg%pq8(kk) * obs%gra(obs%k)
+
+ endif
+
enddo
-
+
+! we apply contribution in grd%variable_ad
+
+ if((drv%chl_assim.eq.1) .or. (drv%multiv.eq.1)) then
+ slicevar(:,1:my_km) = grd%chl_ad(1,:,1:my_km)
+ call ADD_PREVCORE_CONTRIB(ChlExtended_3d, my_km, grd%chl_ad, slicevar(:,1:my_km))
+ ! call ADD_PREVCORE_CONTRIB(ChlExtended_3d, my_km, grd%chl_ad, grd%chl_ad(1,:,:))
+ endif
+ if((bio%N3n.eq.1) .or. (drv%multiv.eq.1)) then
+ slicevar(:,1:grd%km) = grd%n3n_ad(1,:,:)
+ call ADD_PREVCORE_CONTRIB(N3nExtended_3d, grd%km, grd%n3n_ad, slicevar)
+ ! call ADD_PREVCORE_CONTRIB(N3nExtended_3d, grd%km, grd%N3n_ad, grd%n3n_ad(1,:,:))
+ endif
+ if (bio%O2o.eq.1 ) then
+ slicevar(:,1:grd%km) = grd%o2o_ad(1,:,:)
+ call ADD_PREVCORE_CONTRIB(O2oExtended_3d, grd%km, grd%o2o_ad, slicevar)
+ ! call ADD_PREVCORE_CONTRIB(O2oExtended_3d, grd%km, grd%O2o_ad, grd%o2o_ad(1,:,:))
+ endif
+
+
+
end subroutine obs_arg_ad
diff --git a/obs_str.f90 b/obs_str.f90
index 0ab8847..9d5eee7 100644
--- a/obs_str.f90
+++ b/obs_str.f90
@@ -68,6 +68,7 @@ MODULE obs_str
REAL(r8), POINTER :: tim(:) ! Time
REAL(r8), POINTER :: inc(:) ! Increments
REAL(r8), POINTER :: err(:) ! Observational error
+ ! REAL(r8), POINTER :: std(:) ! STD used for EOFs normalization
REAL(r8), POINTER :: res(:) ! residual
INTEGER(i8), POINTER :: ib(:) ! i index of the nearest west point
REAL(r8) , POINTER :: pb(:) ! distance from the nearest west point
diff --git a/obs_vec.f90 b/obs_vec.f90
index 1ca2966..dd9d7fb 100644
--- a/obs_vec.f90
+++ b/obs_vec.f90
@@ -42,7 +42,7 @@ subroutine obs_vec
! -------
! Define observational vector
- obs%no = sat%nc + arg%no
+ obs%no = sat%nc + arg%nc
if(MyId .eq. 0) &
write(drv%dia,*) ' Total number of observations: ', obs%no
@@ -70,7 +70,7 @@ subroutine obs_vec
endif
- ! Observations of chlorophyll
+ ! Observations of satellite chlorophyll
if(drv%sat_obs .eq. 1) then
do i=1,sat%no
if(sat%flc(i).eq.1)then
diff --git a/obsop.f90 b/obsop.f90
index 3a88e4d..466e1f1 100644
--- a/obsop.f90
+++ b/obsop.f90
@@ -42,7 +42,7 @@ subroutine obsop
! ---
! Apply biological repartition of the chlorophyll
- if(drv%chl_assim .eq. 1) &
+ if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) &
call bio_conv
! ---
diff --git a/obsop_ad.f90 b/obsop_ad.f90
index e1d0385..cc90231 100644
--- a/obsop_ad.f90
+++ b/obsop_ad.f90
@@ -52,7 +52,7 @@ subroutine obsop_ad
! ---
! Apply biological repartition of the chlorophyll
- if(drv%chl_assim .eq. 1) &
+ if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) &
call bio_conv_ad
call MPI_Barrier(Var3DCommunicator, ierr)
diff --git a/oceanvar.f90 b/oceanvar.f90
index a635daa..67df41c 100644
--- a/oceanvar.f90
+++ b/oceanvar.f90
@@ -38,6 +38,7 @@ subroutine oceanvar
use set_knd
use drv_str
use mpi_str
+ use da_params
implicit none
@@ -47,7 +48,8 @@ subroutine oceanvar
! ---
! Initialize diagnostics and read namelists
- call def_nml
+ call def_nml ! General DA parameters
+ call def_nml_multi ! DA parameters for multivariate and multiplatform
! ---
! Define grid parameters
@@ -102,10 +104,20 @@ subroutine oceanvar
call wrt_dia
! Write restarts for chl and related variables
- if(drv%chl_assim .eq. 1) &
+ if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) then
call wrt_chl_stat
- if (drv%nut .eq. 1) &
+ ! To write a copy of RSTbefore in RST_after
+ ! In case of assimiation of chl only at some dates
+ if ((drv%nut.eq.0) .and. (NNutVar.gt.0) .and. (drv%multiv.eq.0)) then
+ if (drv%chl_upnut .eq. 0) &
+ call cp_nut_stat
+ if (drv%chl_upnut .eq. 1) &
+ call wrt_upd_nut
+ endif
+ endif
+
+ if ((drv%nut .eq. 1) .or. (drv%multiv.eq.1)) &
call wrt_nut_stat
call sav_itr
diff --git a/rdeofs_chl.f90 b/rdeofs_chl.f90
index ee59cef..9f32a38 100644
--- a/rdeofs_chl.f90
+++ b/rdeofs_chl.f90
@@ -46,7 +46,7 @@ subroutine rdeofs_chl
real(4), allocatable :: x3(:,:,:), x2(:,:)
stat = nf90mpi_open(Var3DCommunicator, trim(EOF_FILE_CHL), NF90_NOWRITE, MPI_INFO_NULL, ncid)
- if (stat /= nf90_noerr) call handle_err("nf90mpi_open", stat)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(EOF_FILE_CHL), stat)
! Get dimensions
stat = nf90mpi_inq_dimid (ncid, 'nreg', idvar)
diff --git a/rdeofs_multi.f90 b/rdeofs_multi.f90
new file mode 100644
index 0000000..2dfcbc7
--- /dev/null
+++ b/rdeofs_multi.f90
@@ -0,0 +1,212 @@
+subroutine rdeofs_multi
+
+ !---------------------------------------------------------------------------
+ ! !
+ ! Copyright 2006 Srdjan Dobricic, CMCC, Bologna !
+ ! !
+ ! This file is part of OceanVar. !
+ ! !
+ ! OceanVar is free software: you can redistribute it and/or modify. !
+ ! it under the terms of the GNU General Public License as published by !
+ ! the Free Software Foundation, either version 3 of the License, or !
+ ! (at your option) any later version. !
+ ! !
+ ! OceanVar is distributed in the hope that it will be useful, !
+ ! but WITHOUT ANY WARRANTY; without even the implied warranty of !
+ ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
+ ! GNU General Public License for more details. !
+ ! !
+ ! You should have received a copy of the GNU General Public License !
+ ! along with OceanVar. If not, see . !
+ ! !
+ !---------------------------------------------------------------------------
+
+ !-----------------------------------------------------------------------
+ ! !
+ ! READ parameters of the MFS_16_72 grid !
+ ! !
+ ! Version 1: S.Dobricic 2006 !
+ ! This routine will have effect only if compiled with netcdf library. !
+ !-----------------------------------------------------------------------
+
+ use set_knd
+ use drv_str
+ use eof_str
+ use grd_str
+ use filenames
+
+ use mpi_str
+ use pnetcdf
+
+ implicit none
+
+ INTEGER(i4) :: stat, ncid, idvar
+ INTEGER(8) :: ik,ireg
+ integer(8) :: neofs, nlevs, nregs, nregsstd
+ integer(KIND=MPI_OFFSET_KIND) :: GlobalStart(3), GlobalCount(3)
+ real(4), allocatable :: std_chl(:),std_n3n(:)
+ real(4), allocatable :: x3(:,:,:), x2(:,:), x1(:)
+
+
+ ! Read std file for normalization
+ stat = nf90mpi_open(Var3DCommunicator, trim(STD_FILE_MULTI), NF90_NOWRITE, MPI_INFO_NULL, ncid)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(STD_FILE_MULTI), stat)
+
+ ! Get dimensions
+ stat = nf90mpi_inq_dimid (ncid, 'nreg', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_dimid nreg", stat)
+ stat = nfmpi_inq_dimlen (ncid, idvar, nregsstd)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_inq_dimlen nregsstd", stat)
+
+ if(MyId .eq. 0) then
+ write(drv%dia,*)'Number of reg in std file : ',nregsstd
+ write(drv%dia,*)'Uses ',ros%neof_multi,' eofs.'
+ endif
+
+ if(ros%nreg .ne. nregsstd) then
+
+ if(MyId .eq. 0) &
+ write(drv%dia,*)'Error: nregsstd differs from nregs'
+
+ call MPI_Abort(Var3DCommunicator, -1, stat)
+
+ endif
+
+ ! Allocate eof arrays and get data
+ ALLOCATE ( std_chl( ros%nreg ) ) ; std_chl = huge(std_chl(1))
+ ALLOCATE ( std_n3n( ros%nreg ) ) ; std_n3n = huge(std_n3n(1))
+ ALLOCATE ( x1( ros%nreg ) )
+ GlobalStart(:) = 1
+ GlobalCount(1) = ros%nreg
+ GlobalCount(2) = ros%kmchl+ros%kmnit
+ GlobalCount(3) = ros%neof_multi
+
+ stat = nf90mpi_inq_varid(ncid, 'stdP_l', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_varid evc", stat)
+ stat = nfmpi_get_vara_real_all(ncid,idvar,GlobalStart(1), GlobalCount(1), x1)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_get_vara_real_all eva", stat)
+
+ std_chl(:) = x1(:)
+
+ stat = nf90mpi_inq_varid(ncid, 'stdN3n', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_varid eva", stat)
+ stat = nfmpi_get_vara_real_all(ncid,idvar,GlobalStart(1), GlobalCount(1), x1)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_get_vara_real_all", stat)
+
+ std_n3n(:) = x1(:)
+
+ ! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST
+ ! ros%evc_multi(:,:,:) = 1.
+ ! ros%eva_multi(:,:) = 1.
+
+ stat = nf90mpi_close(ncid)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_close", stat)
+
+
+
+
+ stat = nf90mpi_open(Var3DCommunicator, trim(EOF_FILE_MULTI), NF90_NOWRITE, MPI_INFO_NULL, ncid)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(EOF_FILE_MULTI), stat)
+
+ ! Get dimensions
+ stat = nf90mpi_inq_dimid (ncid, 'nreg', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_dimid nreg", stat)
+ stat = nfmpi_inq_dimlen (ncid, idvar, nregs)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_inq_dimlen nregs", stat)
+ stat = nf90mpi_inq_dimid (ncid, 'nlev', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_dimid nlev", stat)
+ stat = nfmpi_inq_dimlen (ncid, idvar, nlevs)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_inq_dimlen nlevs", stat)
+ stat = nf90mpi_inq_dimid (ncid, 'neof', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_dimid neof", stat)
+ stat = nfmpi_inq_dimlen (ncid, idvar, len = neofs)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_inq_dimlen neofs", stat)
+
+ if(MyId .eq. 0) then
+ write(drv%dia,*)'Eof dimensions for multivariate are: ',ros%nreg, ros%kmchl+ros%kmnit, neofs
+ write(drv%dia,*)'Uses ',ros%neof_multi,' eofs.'
+ endif
+
+ if(ros%nreg .ne. nregs) then
+
+ if(MyId .eq. 0) &
+ write(drv%dia,*)'Error: ros%nreg differs from nregs'
+
+ call MPI_Abort(Var3DCommunicator, -1, stat)
+
+ endif
+
+ if(ros%neof_multi .gt. neofs) then
+
+ if(MyId .eq. 0) &
+ write(drv%dia,*)'Error: Requires more Eofs than available in the input file.'
+ call MPI_Abort(Var3DCommunicator, -1, stat)
+
+ else if(ros%neof_multi .lt. neofs) then
+
+ if(MyId .eq. 0) then
+ write(drv%dia,*)'Warning: ros%neof_multi < neofs!'
+ write(drv%dia,*)'ros%neof_multi =', ros%neof_multi
+ write(drv%dia,*)'neofs =', neofs
+ write(drv%dia,*)'continue using ros%neof_multi'
+ write(*,*)'Warning: ros%neof_multi < neofs!'
+ write(*,*)'ros%neof_multi =', ros%neof_multi
+ write(*,*)'neofs =', neofs
+ write(*,*)'continue using ros%neof_multi'
+ endif
+ endif
+
+ if((ros%kmchl+ros%kmnit) .ne. nlevs) then
+ if(MyId .eq. 0) &
+ write(drv%dia,*)'Error: Vertical dimension different than in the input file.'
+
+ call MPI_Abort(Var3DCommunicator, -1, stat)
+ endif
+
+ ! Allocate eof arrays and get data
+ ALLOCATE ( ros%evc_multi( ros%nreg, ros%kmchl+ros%kmnit, ros%neof_multi) ) ; ros%evc_multi = huge(ros%evc_multi(1,1,1))
+ ALLOCATE ( ros%eva_multi( ros%nreg, ros%neof_multi) ) ; ros%eva_multi = huge(ros%eva_multi(1,1))
+ ALLOCATE ( x3( ros%nreg, ros%kmchl+ros%kmnit, ros%neof_multi) )
+ ALLOCATE ( x2( ros%nreg, ros%neof_multi) )
+ GlobalStart(:) = 1
+ GlobalCount(1) = ros%nreg
+ GlobalCount(2) = ros%kmchl+ros%kmnit
+ GlobalCount(3) = ros%neof_multi
+
+ stat = nf90mpi_inq_varid(ncid, 'evc', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_varid evc", stat)
+ stat = nfmpi_get_vara_real_all(ncid,idvar,GlobalStart, GlobalCount, x3)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_get_vara_real_all eva", stat)
+
+ ros%evc_multi(:,:,:) = x3(:,:,:)
+
+ do ireg=1,ros%nreg
+ do ik=1,ros%kmchl
+ ros%evc_multi(ireg,ik,:) = ros%evc_multi(ireg,ik,:)*std_chl(ireg)
+ enddo
+ do ik=ros%kmchl+1,ros%kmchl+ros%kmnit
+ ros%evc_multi(ireg,ik,:) = ros%evc_multi(ireg,ik,:)*std_n3n(ireg)
+ enddo
+ enddo
+
+ GlobalCount(1) = ros%nreg
+ GlobalCount(2) = ros%neof_multi
+
+ stat = nf90mpi_inq_varid(ncid, 'eva', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_varid eva", stat)
+ stat = nfmpi_get_vara_real_all(ncid,idvar,GlobalStart(1:2), GlobalCount(1:2), x2)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_get_vara_real_all", stat)
+ ros%eva_multi(:,:) = x2(:,:)
+
+ ! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST
+ ! ros%evc_multi(:,:,:) = 1.
+ ! ros%eva_multi(:,:) = 1.
+
+ stat = nf90mpi_close(ncid)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_close", stat)
+
+ DEALLOCATE(x3, x2, std_chl, std_n3n)
+
+end subroutine rdeofs_multi
+
+
diff --git a/rdeofs_n3n.f90 b/rdeofs_n3n.f90
index a45175c..5e158f2 100644
--- a/rdeofs_n3n.f90
+++ b/rdeofs_n3n.f90
@@ -47,7 +47,7 @@ subroutine rdeofs_n3n
stat = nf90mpi_open(Var3DCommunicator, trim(EOF_FILE_N3N), NF90_NOWRITE, MPI_INFO_NULL, ncid)
- if (stat /= nf90_noerr) call handle_err("nf90mpi_open", stat)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(EOF_FILE_N3N), stat)
! Get dimensions
stat = nf90mpi_inq_dimid (ncid, 'nreg', idvar)
diff --git a/rdeofs_o2o.f90 b/rdeofs_o2o.f90
index c6c840f..83290b8 100644
--- a/rdeofs_o2o.f90
+++ b/rdeofs_o2o.f90
@@ -46,7 +46,7 @@ subroutine rdeofs_o2o
real(4), allocatable :: x3(:,:,:), x2(:,:)
stat = nf90mpi_open(Var3DCommunicator, trim(EOF_FILE_O2O), NF90_NOWRITE, MPI_INFO_NULL, ncid)
- if (stat /= nf90_noerr) call handle_err("nf90mpi_open", stat)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(EOF_FILE_O2O), stat)
! Get dimensions
stat = nf90mpi_inq_dimid (ncid, 'nreg', idvar)
diff --git a/rdrcorr.f90 b/rdrcorr.f90
index cbddd8e..ca4a485 100644
--- a/rdrcorr.f90
+++ b/rdrcorr.f90
@@ -38,9 +38,9 @@ subroutine rdrcorr
integer(KIND=MPI_OFFSET_KIND) :: GlobalStart(3), GlobalCount(3)
real(r4), ALLOCATABLE :: x3(:,:,:)
- !write(*,*)trim(RCORR_FILE)
+ ! print*,RCORR_FILE
stat = nf90mpi_open(Var3DCommunicator, trim(RCORR_FILE), NF90_NOWRITE, MPI_INFO_NULL, ncid)
- if (stat /= nf90_noerr) call handle_err("nf90mpi_open",stat)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(RCORR_FILE),stat)
ALLOCATE ( x3(GlobalRow,GlobalCol,grd%km))
GlobalStart(:) = 1
diff --git a/readAnisotropy.f90 b/readAnisotropy.f90
index 4a4b4cc..db9d3f3 100644
--- a/readAnisotropy.f90
+++ b/readAnisotropy.f90
@@ -38,7 +38,7 @@ subroutine readAnisotropy
real(r4), ALLOCATABLE :: x2(:,:)
stat = nf90mpi_open(Var3DCommunicator, trim(ANIS_FILE), NF90_NOWRITE, MPI_INFO_NULL, ncid)
- if (stat /= nf90_noerr) call handle_err("nf90mpi_open",stat)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(ANIS_FILE),stat)
ALLOCATE ( x2(GlobalRow,GlobalCol))
GlobalStart(:) = 1
@@ -83,4 +83,4 @@ subroutine readAnisotropy
if (stat /= nf90_noerr) call handle_err("nf90mpi_close", stat)
DEALLOCATE(x2)
-end subroutine readAnisotropy
\ No newline at end of file
+end subroutine readAnisotropy
diff --git a/readChlNutCov.f90 b/readChlNutCov.f90
new file mode 100644
index 0000000..b39898b
--- /dev/null
+++ b/readChlNutCov.f90
@@ -0,0 +1,127 @@
+subroutine readChlNutCov
+
+!---------------------------------------------------------------------------
+! !
+! Copyright 2015 Anna teruzzi, OGS Trieste !
+! !
+! This file is part of OceanVar. !
+! !
+! OceanVar is free software: you can redistribute it and/or modify. !
+! it under the terms of the GNU General Public License as published by !
+! the Free Software Foundation, either version 3 of the License, or !
+! (at your option) any later version. !
+! !
+! OceanVar is distributed in the hope that it will be useful, !
+! but WITHOUT ANY WARRANTY; without even the implied warranty of !
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
+! GNU General Public License for more details. !
+! !
+! You should have received a copy of the GNU General Public License !
+! along with OceanVar. If not, see . !
+! !
+!---------------------------------------------------------------------------
+
+
+!---------------------------------------------------------------------------
+! read covariances between Nitrate and Phsosphate !
+! to update Phosphate with assimilation of Nitrate !
+!---------------------------------------------------------------------------
+
+
+ use set_knd
+ use drv_str
+ use grd_str
+ use bio_str
+ use filenames
+
+ use mpi_str
+ use pnetcdf
+
+ implicit none
+
+ integer(i4) :: stat, ncid, idvar
+ integer(i4) :: i,j,k
+ !integer(KIND=MPI_OFFSET_KIND) :: GlobalStart(3), GlobalCount(3)
+ real(r4), ALLOCATABLE :: x3(:,:,:)
+
+
+ !write(*,*)trim(RCORR_FILE)
+ stat = nf90mpi_open(Var3DCommunicator, trim(NUTCHLCOV_FILE), NF90_NOWRITE, MPI_INFO_NULL, ncid)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(NUTCHLCOV_FILE),stat)
+
+ ALLOCATE ( x3(grd%im, grd%jm, grd%km))
+ ALLOCATE ( bio%covn3n_chl(grd%im, grd%jm, grd%km)); bio%covn3n_chl(:,:,:) = 0.0
+ ALLOCATE ( bio%covn1p_chl(grd%im, grd%jm, grd%km)); bio%covn1p_chl(:,:,:) = 0.0
+
+ ! covn3n_chl
+ x3(:,:,:) = 0.0
+
+ stat = nf90mpi_inq_varid (ncid, 'covp_l_n3n', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_varid radius",stat)
+ stat = nfmpi_get_vara_real_all (ncid, idvar, MyStart, MyCount, x3)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_get_vara_real_all radius",stat)
+
+ do k=1,grd%km
+ do j=1,grd%jm
+ do i=1,grd%im
+ if(x3(i,j,k) .lt. 1.e20) then
+
+ bio%covn3n_chl(i,j,k) = bio%covn3n_chl(i,j,k) + x3(i,j,k)
+
+ else
+ bio%covn3n_chl(i,j,k) = x3(i,j,k)
+ if(grd%msk(i,j,k) .eq. 1) then
+ write(*,*) "Warning!! Bad mask point in N3n N1p covaraince!"
+ write(*,*) "i=",i," j=",j," k=",k
+ write(*,*) "grd%msk(i,j,k)=",grd%msk
+ write(*,*) "bio%covn3n_chl(i,j,k)=",bio%covn3n_chl(i,j,k)
+ write(*,*) "Aborting.."
+ call MPI_Abort(Var3DCommunicator, -1, stat)
+ endif
+
+ endif
+
+ enddo
+ enddo
+ enddo
+
+
+ ! covn1p_chl
+ x3(:,:,:) = 0.0
+
+ stat = nf90mpi_inq_varid (ncid, 'covp_l_n1p', idvar)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_varid radius",stat)
+ stat = nfmpi_get_vara_real_all (ncid, idvar, MyStart, MyCount, x3)
+ if (stat /= nf90_noerr) call handle_err("nfmpi_get_vara_real_all radius",stat)
+
+ do k=1,grd%km
+ do j=1,grd%jm
+ do i=1,grd%im
+ if(x3(i,j,k) .lt. 1.e20) then
+
+ bio%covn1p_chl(i,j,k) = bio%covn1p_chl(i,j,k) + x3(i,j,k)
+
+ else
+ bio%covn1p_chl(i,j,k) = x3(i,j,k)
+ if(grd%msk(i,j,k) .eq. 1) then
+ write(*,*) "Warning!! Bad mask point in N3n N1p covaraince!"
+ write(*,*) "i=",i," j=",j," k=",k
+ write(*,*) "grd%msk(i,j,k)=",grd%msk
+ write(*,*) "bio%covn1p_chl(i,j,k)=",bio%covn1p_chl(i,j,k)
+ write(*,*) "Aborting.."
+ call MPI_Abort(Var3DCommunicator, -1, stat)
+ endif
+
+ endif
+
+ enddo
+ enddo
+ enddo
+
+
+ stat = nf90mpi_close(ncid)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_close", stat)
+
+ DEALLOCATE(x3)
+
+end subroutine readChlNutCov
diff --git a/readChlStat.f90 b/readChlStat.f90
index 5bec730..4e0437b 100644
--- a/readChlStat.f90
+++ b/readChlStat.f90
@@ -42,7 +42,7 @@ subroutine readChlStat
INTEGER(i4) :: ncid, VarId, ierr, iVar
INTEGER(i4) :: i,j,k,l,m
- CHARACTER(LEN=51) :: RstFileName
+ CHARACTER(LEN=47) :: RstFileName
CHARACTER(LEN=3) :: MyVarName
REAL(4), ALLOCATABLE :: x3(:,:,:)
@@ -60,7 +60,7 @@ subroutine readChlStat
if(iVar .gt. NPhytoVar) cycle
MyVarName = DA_VarList(iVar)
- RstFileName = 'DA__FREQ_1/RST.'//ShortDate//'.'//MyVarName//'.nc'
+ RstFileName = 'DA__FREQ_1/RSTbefore.'//ShortDate//'.'//MyVarName//'.nc'
if(drv%Verbose .eq. 1) then
if(MyId .eq. 0) &
@@ -92,7 +92,7 @@ subroutine readChlStat
if(grd%msk(i,j,k) .eq. 1) then
write(*,*) "Warning!! Bad mask point in bio structure!"
write(*,*) "i=",i," j=",j," k=",k
- write(*,*) "grd%msk(i,j,k)=",grd%msk
+ write(*,*) "grd%msk(i,j,k)=",grd%msk(i,j,k)
write(*,*) "bio%InitialChl(i,j,k)=",bio%InitialChl(i,j,k)
write(*,*) "Aborting.."
call MPI_Abort(Var3DCommunicator, -1, ierr)
diff --git a/readGrid.f90 b/readGrid.f90
index 9761768..992388e 100644
--- a/readGrid.f90
+++ b/readGrid.f90
@@ -3,6 +3,7 @@ subroutine readGrid
use set_knd
use drv_str
use grd_str
+ use eof_str
use filenames
use mpi_str
use pnetcdf
@@ -12,7 +13,7 @@ subroutine readGrid
implicit none
- integer(i8) :: ierr, ncid
+ integer(i8) :: ierr, ncid, my_km
integer(i8) :: jpreci, jprecj
integer(i8) :: VarId
real(r4), ALLOCATABLE :: x3(:,:,:), x2(:,:), x1(:)
@@ -101,6 +102,20 @@ subroutine readGrid
ALLOCATE(SendTop(grd%jm), RecBottom(grd%jm))
ALLOCATE(SendBottom(grd%jm), RecTop(grd%jm))
+ ALLOCATE(SendTop_2d( grd%jm,grd%km), SendBottom_2d (grd%jm,grd%km))
+ ALLOCATE(RecBottom_2d( grd%jm,grd%km), RecTop_2d( grd%jm,grd%km))
+ if(drv%multiv.eq.0) then
+ ALLOCATE(ChlExtended_3d (grd%im+1, grd%jm, grd%km))
+ else if(drv%multiv.eq.1) then
+ ALLOCATE(ChlExtended_3d (grd%im+1, grd%jm, ros%kmchl))
+ end if
+ ALLOCATE(N3nExtended_3d (grd%im+1, grd%jm, grd%km))
+ ALLOCATE(O2oExtended_3d (grd%im+1, grd%jm, grd%km))
+
+
+
+
+
ALLOCATE ( grd%reg(grd%im,grd%jm)) ; grd%reg = huge(grd%reg(1,1))
ALLOCATE ( grd%msk(grd%im,grd%jm,grd%km)) ; grd%msk = huge(grd%msk(1,1,1))
ALLOCATE ( grd%dep(grd%km)) ; grd%dep = huge(grd%dep(1))
@@ -120,23 +135,33 @@ subroutine readGrid
ALLOCATE ( grd%inx(GlobalRow,localCol,grd%km)) ; grd%inx = huge(grd%inx(1,1,1))
ALLOCATE ( grd%jnx(localRow,GlobalCol,grd%km)) ; grd%jnx = huge(grd%jnx(1,1,1))
- if(drv%chl_assim .eq. 1) then
- ALLOCATE ( grd%chl(grd%im,grd%jm,grd%km) ) ; grd%chl = huge(grd%chl(1,1,1))
- ALLOCATE ( grd%chl_ad(grd%im,grd%jm,grd%km) ) ; grd%chl_ad = huge(grd%chl_ad(1,1,1))
- endif
- if(drv%nut .eq. 1) then
- if(bio%N3n .eq. 1) then
- ALLOCATE ( grd%n3n(grd%im,grd%jm,grd%km) ) ; grd%n3n = huge(grd%n3n(1,1,1))
- ALLOCATE ( grd%n3n_ad(grd%im,grd%jm,grd%km) ) ; grd%n3n_ad = huge(grd%n3n_ad(1,1,1))
+ if(drv%multiv .eq. 0) then
+ if(drv%chl_assim .eq. 1) then
+ ALLOCATE ( grd%chl(grd%im,grd%jm,grd%km) ) ; grd%chl = huge(grd%chl(1,1,1))
+ ALLOCATE ( grd%chl_ad(grd%im,grd%jm,grd%km) ) ; grd%chl_ad = huge(grd%chl_ad(1,1,1))
+ ALLOCATE ( bio%phy(grd%im,grd%jm,grd%km,bio%nphy,bio%ncmp) ) ; bio%phy = huge(bio%phy(1,1,1,1,1))
+ ALLOCATE ( bio%phy_ad(grd%im,grd%jm,grd%km,bio%nphy,bio%ncmp) ) ; bio%phy_ad = huge(bio%phy_ad(1,1,1,1,1))
endif
- if(bio%O2o .eq. 1) then
- ALLOCATE ( grd%o2o(grd%im,grd%jm,grd%km) ) ; grd%o2o = huge(grd%o2o(1,1,1))
- ALLOCATE ( grd%o2o_ad(grd%im,grd%jm,grd%km) ) ; grd%o2o_ad = huge(grd%o2o_ad(1,1,1))
+ if(drv%nut .eq. 1) then
+ if(bio%N3n .eq. 1) then
+ ALLOCATE ( grd%n3n(grd%im,grd%jm,grd%km) ) ; grd%n3n = huge(grd%n3n(1,1,1))
+ ALLOCATE ( grd%n3n_ad(grd%im,grd%jm,grd%km) ) ; grd%n3n_ad = huge(grd%n3n_ad(1,1,1))
+ endif
+ if(bio%O2o .eq. 1) then
+ ALLOCATE ( grd%o2o(grd%im,grd%jm,grd%km) ) ; grd%o2o = huge(grd%o2o(1,1,1))
+ ALLOCATE ( grd%o2o_ad(grd%im,grd%jm,grd%km) ) ; grd%o2o_ad = huge(grd%o2o_ad(1,1,1))
+ endif
endif
+
+ else if(drv%multiv .eq. 1) then
+ ALLOCATE ( grd%chl(grd%im,grd%jm,ros%kmchl) ) ; grd%chl = huge(grd%chl(1,1,1))
+ ALLOCATE ( grd%chl_ad(grd%im,grd%jm,ros%kmchl) ) ; grd%chl_ad = huge(grd%chl_ad(1,1,1))
+ ALLOCATE ( grd%n3n(grd%im,grd%jm,ros%kmnit) ) ; grd%n3n = huge(grd%n3n(1,1,1))
+ ALLOCATE ( grd%n3n_ad(grd%im,grd%jm,ros%kmnit) ) ; grd%n3n_ad = huge(grd%n3n_ad(1,1,1))
+ ALLOCATE ( bio%phy(grd%im,grd%jm,ros%kmchl,bio%nphy,bio%ncmp) ) ; bio%phy = huge(bio%phy(1,1,1,1,1))
+ ALLOCATE ( bio%phy_ad(grd%im,grd%jm,ros%kmchl,bio%nphy,bio%ncmp) ) ; bio%phy_ad = huge(bio%phy_ad(1,1,1,1,1))
endif
- ALLOCATE ( bio%phy(grd%im,grd%jm,grd%km,bio%nphy,bio%ncmp) ) ; bio%phy = huge(bio%phy(1,1,1,1,1))
- ALLOCATE ( bio%phy_ad(grd%im,grd%jm,grd%km,bio%nphy,bio%ncmp) ) ; bio%phy_ad = huge(bio%phy_ad(1,1,1,1,1))
ALLOCATE ( x3(grd%im,grd%jm,grd%km)) ; x3 = huge(x3(1,1,1))
ALLOCATE ( x2(grd%im,grd%jm)) ; x2 = huge(x2(1,1))
@@ -208,7 +233,6 @@ subroutine readGrid
! *****************************************************************************************
! *****************************************************************************************
- call CreateMpiWindows
end subroutine readGrid
@@ -217,13 +241,14 @@ subroutine DomainDecomposition
use drv_str
use mpi_str
use grd_str
+ use eof_str
implicit none
integer, allocatable :: ilcit(:,:), ilcjt(:,:), BalancedSlice(:,:)
integer(i8) :: ji, jj, TmpInt, ierr ! jpi, jpj, nn, i
integer(i8) :: GlobalRestCol, GlobalRestRow
- integer(i8) :: i, j, k, kk
+ integer(i8) :: i, j, k, kk, my_km
integer(i8) :: NCoastX, NCoastY, TmpCoast
integer(i8) :: NRows, NCols
integer(i8) :: SliceRestRow, SliceRestCol
@@ -312,10 +337,16 @@ subroutine DomainDecomposition
SendDisplX3D(1) = 0
RecDisplX3D(1) = 0
+ SendDisplX3D_chl(1) = 0
+ RecDisplX3D_chl(1) = 0
SendDisplX2D(1) = 0
RecDisplX2D(1) = 0
+ my_km = grd%km
+ if(drv%multiv.eq.1) &
+ my_km = ros%kmchl
+
do i=1,NumProcI
if(i-1 .lt. SliceRestCol) then
OffsetRow = 1
@@ -332,6 +363,9 @@ subroutine DomainDecomposition
SendCountX3D(i) = (grd%jm / NumProcI + OffsetRow) * grd%im * grd%km
RecCountX3D(i) = localCol * grd%km * (GlobalRow / NumProcI + OffsetCol)
+ SendCountX3D_chl(i) = (grd%jm / NumProcI + OffsetRow) * grd%im * my_km
+ RecCountX3D_chl(i) = localCol * my_km * (GlobalRow / NumProcI + OffsetCol)
+
SendCountX2D(i) = (grd%jm / NumProcI + OffsetRow) * grd%im
RecCountX2D(i) = localCol * (GlobalRow / NumProcI + OffsetCol)
@@ -339,6 +373,9 @@ subroutine DomainDecomposition
SendDisplX3D(i+1) = SendDisplX3D(i) + SendCountX3D(i)
RecDisplX3D(i+1) = RecDisplX3D(i) + RecCountX3D(i)
+ SendDisplX3D_chl(i+1) = SendDisplX3D_chl(i) + SendCountX3D_chl(i)
+ RecDisplX3D_chl(i+1) = RecDisplX3D_chl(i) + RecCountX3D_chl(i)
+
SendDisplX2D(i+1) = SendDisplX2D(i) + SendCountX2D(i)
RecDisplX2D(i+1) = RecDisplX2D(i) + RecCountX2D(i)
end if
@@ -441,37 +478,4 @@ subroutine handle_err(err_msg, errcode)
end subroutine handle_err
-subroutine CreateMpiWindows
-
- use grd_str
- use mpi_str
- use drv_str
- use bio_str
-
- implicit none
- ! include "mpif.h"
-
- integer :: ierr
- integer(kind=MPI_ADDRESS_KIND) :: nbytes, lenreal, dummy
-
- ! lenreal = 8
- call MPI_Type_get_extent(MPI_REAL8, dummy, lenreal, ierr)
- nbytes = grd%im*grd%jm*grd%km*lenreal
-
- if(drv%chl_assim .eq. 1) then
- call MPI_Win_create(grd%chl, nbytes, lenreal, MPI_INFO_NULL, Var3DCommunicator, MpiWinChl, ierr)
- call MPI_Win_create(grd%chl_ad, nbytes, lenreal, MPI_INFO_NULL, Var3DCommunicator, MpiWinChlAd, ierr)
- endif
- if(drv%nut .eq. 1) then
- if(bio%N3n .eq. 1) then
- call MPI_Win_create(grd%n3n, nbytes, lenreal, MPI_INFO_NULL, Var3DCommunicator, MpiWinN3n, ierr)
- call MPI_Win_create(grd%n3n_ad, nbytes, lenreal, MPI_INFO_NULL, Var3DCommunicator, MpiWinN3nAd, ierr)
- endif
- if(bio%O2o .eq. 1) then
- call MPI_Win_create(grd%o2o, nbytes, lenreal, MPI_INFO_NULL, Var3DCommunicator, MpiWinO2o, ierr)
- call MPI_Win_create(grd%o2o_ad, nbytes, lenreal, MPI_INFO_NULL, Var3DCommunicator, MpiWinO2oAd, ierr)
- endif
- endif
-
-end subroutine CreateMpiWindows
diff --git a/readNutCov.f90 b/readNutCov.f90
index 2f1427d..1595d71 100644
--- a/readNutCov.f90
+++ b/readNutCov.f90
@@ -44,16 +44,15 @@ subroutine readNutCov
!integer(KIND=MPI_OFFSET_KIND) :: GlobalStart(3), GlobalCount(3)
real(r4), ALLOCATABLE :: x3(:,:,:)
-
- !write(*,*)trim(RCORR_FILE)
- stat = nf90mpi_open(Var3DCommunicator, trim(NUTCOV_FILE), NF90_NOWRITE, MPI_INFO_NULL, ncid)
- if (stat /= nf90_noerr) call handle_err("nf90mpi_open",stat)
-
ALLOCATE ( x3(grd%im, grd%jm, grd%km))
ALLOCATE ( bio%covn3n_n1p(grd%im, grd%jm, grd%km)); bio%covn3n_n1p(:,:,:) = 0.0
x3(:,:,:) = 0.0
+ !write(*,*)trim(RCORR_FILE)
+ stat = nf90mpi_open(Var3DCommunicator, trim(NUTCOV_FILE), NF90_NOWRITE, MPI_INFO_NULL, ncid)
+ if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(NUTCOV_FILE),stat)
+
stat = nf90mpi_inq_varid (ncid, 'covn3n_n1p', idvar)
if (stat /= nf90_noerr) call handle_err("nf90mpi_inq_varid radius",stat)
stat = nfmpi_get_vara_real_all (ncid, idvar, MyStart, MyCount, x3)
diff --git a/readNutStat.f90 b/readNutStat.f90
index 050d640..b035c8c 100644
--- a/readNutStat.f90
+++ b/readNutStat.f90
@@ -42,7 +42,7 @@ subroutine readNutStat
INTEGER(i4) :: ncid, VarId, ierr, iVar
INTEGER(i4) :: i,j,k,l
- CHARACTER(LEN=51) :: RstFileName
+ CHARACTER(LEN=47) :: RstFileName
CHARACTER(LEN=3) :: MyVarName
REAL(4), ALLOCATABLE :: x3(:,:,:)
@@ -62,7 +62,7 @@ subroutine readNutStat
endif
MyVarName = DA_VarList(iVar)
- RstFileName = 'DA__FREQ_1/RST.'//ShortDate//'.'//MyVarName//'.nc'
+ RstFileName = 'DA__FREQ_1/RSTbefore.'//ShortDate//'.'//MyVarName//'.nc'
if(drv%Verbose .eq. 1) then
if(MyId .eq. 0) &
@@ -90,8 +90,8 @@ subroutine readNutStat
if(grd%msk(i,j,k) .eq. 1) then
write(*,*) "Warning!! Bad mask point in bio structure!"
write(*,*) "i=",i," j=",j," k=",k
- write(*,*) "grd%msk(i,j,k)=",grd%msk
- write(*,*) "bio%InitialChl(i,j,k)=",bio%InitialChl(i,j,k)
+ write(*,*) "grd%msk(i,j,k)=",grd%msk(i,j,k)
+ write(*,*) "bio%InitialNut(i,j,k)=",bio%InitialNut(i,j,k,l)
write(*,*) "Aborting.."
call MPI_Abort(Var3DCommunicator, -1, ierr)
endif
diff --git a/res_inc.f90 b/res_inc.f90
index e50a492..edc891b 100644
--- a/res_inc.f90
+++ b/res_inc.f90
@@ -33,10 +33,26 @@ subroutine res_inc
use drv_str
use grd_str
use obs_str
+ use bio_str
implicit none
- grd%chl_ad(:,:,:) = 0.0 ! OMP
+ if (drv%multiv .eq. 0) then
+ if (drv%chl_assim .eq. 1) then
+ grd%chl_ad(:,:,:) = 0.0 ! OMP
+ end if
+
+ if (drv%nut .eq. 1) then
+ if (bio%n3n .eq. 1) &
+ grd%n3n_ad(:,:,:) = 0.0
+ if (bio%o2o .eq. 1) &
+ grd%o2o_ad(:,:,:) = 0.0
+ endif
+
+ else if(drv%multiv .eq.1) then
+ grd%chl_ad(:,:,:) = 0.0 ! OMP
+ grd%n3n_ad(:,:,:) = 0.0
+ endif
obs%gra(:) = obs%amo(:) / obs%err(:) ! OMP
diff --git a/resid.f90 b/resid.f90
index 9208be9..faba056 100644
--- a/resid.f90
+++ b/resid.f90
@@ -50,9 +50,10 @@ subroutine resid
endif
enddo
endif
+
! ---
- ! Observations of chlorophyll
+ ! Observations of satellite chlorophyll
if(drv%sat_obs .eq. 1) then
do i=1,sat%no
if(sat%flc(i).eq.1)then
diff --git a/sav_itr.f90 b/sav_itr.f90
index 078fc6c..37e6398 100644
--- a/sav_itr.f90
+++ b/sav_itr.f90
@@ -40,11 +40,10 @@ subroutine sav_itr
use rcfl
use mpi_str
use bio_str
+ use da_params
implicit none
- ! free MPI RMA Windows
- call FreeWindows
! ---
! Grid structure
@@ -67,6 +66,8 @@ subroutine sav_itr
DEALLOCATE( grd%bey)
! Chlorophyll vectors
+ if(drv%multiv.eq.0) then
+
if(drv%chl_assim .eq. 1) then
DEALLOCATE( grd%chl)
DEALLOCATE( grd%chl_ad)
@@ -82,6 +83,15 @@ subroutine sav_itr
endif
endif
+ endif
+
+ if(drv%multiv.eq.1) then
+ DEALLOCATE( grd%chl)
+ DEALLOCATE( grd%chl_ad)
+ DEALLOCATE( grd%n3n)
+ DEALLOCATE( grd%n3n_ad)
+ endif
+
! Observational vector
DEALLOCATE( obs%inc, obs%amo, obs%res)
DEALLOCATE( obs%err, obs%gra)
@@ -105,14 +115,32 @@ subroutine sav_itr
DEALLOCATE( ctl%x_c, ctl%g_c)
! Bio structure
- if(drv%chl_assim .eq. 1) then
+ if(drv%multiv.eq.0) then
+ if(drv%chl_assim .eq. 1) then
+ DEALLOCATE( bio%phy, bio%phy_ad)
+ DEALLOCATE( bio%cquot, bio%pquot)
+ DEALLOCATE( bio%InitialChl)
+ if ((drv%nut .eq. 0) .and. (NNutVar .gt. 0)) then
+ DEALLOCATE( bio%InitialNut)
+ if (drv%chl_upnut.eq.1) then
+ ! DEALLOCATE( bio%covn3n_n1p)
+ DEALLOCATE( bio%covn1p_chl)
+ DEALLOCATE( bio%covn3n_chl)
+ endif
+ endif
+ endif
+ if(drv%nut .eq. 1) then
+ DEALLOCATE( bio%InitialNut)
+ if(bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) DEALLOCATE( bio%covn3n_n1p)
+ endif
+ endif
+
+ if(drv%multiv.eq.1) then
DEALLOCATE( bio%phy, bio%phy_ad)
DEALLOCATE( bio%cquot, bio%pquot)
DEALLOCATE( bio%InitialChl)
- endif
- if(drv%nut .eq. 1) then
DEALLOCATE( bio%InitialNut)
- DEALLOCATE( bio%covn3n_n1p)
+ if(bio%updateN1p.eq.1) DEALLOCATE( bio%covn3n_n1p)
endif
DEALLOCATE(SurfaceWaterPoints)
@@ -127,35 +155,10 @@ subroutine sav_itr
DEALLOCATE( bta_rcx)
DEALLOCATE( alp_rcy)
DEALLOCATE( bta_rcy)
+ DEALLOCATE( grd%global_msk)
if(MyId .eq. 0) write(*,*) ' DEALLOCATION DONE'
end subroutine sav_itr
-subroutine FreeWindows
-
- use grd_str
- use mpi_str
- use drv_str
- use bio_str
-
- implicit none
-
- integer ierr
-
- if(drv%chl_assim .eq. 1) then
- call MPI_Win_free(MpiWinChl, ierr)
- call MPI_Win_free(MpiWinChlAd, ierr)
- endif
- if(drv%nut .eq. 1) then
- if(bio%n3n .eq. 1) then
- call MPI_Win_free(MpiWinN3n, ierr)
- call MPI_Win_free(MpiWinN3nAd, ierr)
- endif
- if(bio%o2o .eq. 1) then
- call MPI_Win_free(MpiWinO2o, ierr)
- call MPI_Win_free(MpiWinO2oAd, ierr)
- endif
- endif
-end subroutine FreeWindows
\ No newline at end of file
diff --git a/tao_minimizer.f90 b/tao_minimizer.f90
index 93204f6..44fcee6 100644
--- a/tao_minimizer.f90
+++ b/tao_minimizer.f90
@@ -1,6 +1,11 @@
subroutine tao_minimizer
+#include "petscversion.h"
+#if PETSC_VERSION_GE(3,8,0)
+#include "petsc/finclude/petscvec.h"
+#else
#include "petsc/finclude/petscvecdef.h"
+#endif
use drv_str
use ctl_str
@@ -14,7 +19,7 @@ subroutine tao_minimizer
PetscErrorCode :: ierr
Tao :: tao
Vec :: MyState ! array that stores the (temporary) state
- PetscInt :: n, M, GlobalStart, MyEnd, iter
+ PetscInt :: n, M, GlobalStart, MyEnd, iter!, maxfeval
PetscReal :: fval, gnorm, cnorm, xdiff
PetscScalar :: MyTolerance
TaoConvergedReason :: reason
@@ -89,7 +94,12 @@ subroutine tao_minimizer
! Set initial solution array, MyBounds and MyFuncAndGradient routines
call TaoSetInitialVector(tao, MyState, ierr)
CHKERRQ(ierr)
+#include
+#if PETSC_VERSION_GE(3,8,0)
+ call TaoSetObjectiveAndGradientRoutine(tao, MyFuncAndGradient, PETSC_NULL_VEC, ierr)
+#else
call TaoSetObjectiveAndGradientRoutine(tao, MyFuncAndGradient, PETSC_NULL_OBJECT, ierr)
+#endif
CHKERRQ(ierr)
! Calling costf in order to compute
@@ -108,14 +118,21 @@ subroutine tao_minimizer
if(MyId .eq. 0) then
print*, "Setting MyTolerance", MyTolerance
+ print*, "------- MaxGrad", MaxGrad
print*, ""
write(drv%dia,*) "Setting MyTolerance", MyTolerance
+ write(drv%dia,*) "------- MaxGrad", MaxGrad
endif
! setting tolerances
call TaoSetTolerances(tao, MyTolerance, 1.d-4, ctl%pgper, ierr)
CHKERRQ(ierr)
+ ! setting max number of fucntion evaluation
+ !maxfeval = 300
+ call TaoSetMaximumFunctionEvaluations(tao, 100, ierr)
+ CHKERRQ(ierr)
+
! calling the solver to minimize the problem
call TaoSolve(tao, ierr)
CHKERRQ(ierr)
@@ -131,13 +148,33 @@ subroutine tao_minimizer
call TaoGetSolutionStatus(tao, iter, fval, gnorm, cnorm, xdiff, reason, ierr)
if(reason .lt. 0) then
- if(MyId .eq. 0) then
- print*, "TAO failed to find a solution"
- print*, "Aborting.."
- endif
- call MPI_Barrier(Var3DCommunicator, ierr)
- call MPI_Abort(Var3DCommunicator, -1, ierr)
- endif
+
+ if( ((reason.eq.-6) .or. (reason.eq.-5)) .and. (drv%MyCounter .gt. 12) ) then
+ if(MyId .eq. 0) then
+ print*, "TAO failed to find a solution"
+ print*, "fval..", fval
+ print*, "gnorm.", gnorm
+ print*, "reason", reason
+ print*, "iter", iter
+ print*, " MyCount", drv%MyCounter
+ print*, "BUT assigning a solution "
+ endif
+
+ else
+ if(MyId .eq. 0) then
+ print*, "TAO failed to find a solution"
+ print*, "fval..", fval
+ print*, "gnorm.", gnorm
+ print*, "reason", reason
+ print*, "iter", iter
+ print*, " MyCount", drv%MyCounter
+ print*, "Aborting.."
+ endif
+ call MPI_Barrier(Var3DCommunicator, ierr)
+ call MPI_Abort(Var3DCommunicator, -1, ierr)
+ endif ! reason -6 or -5
+
+ endif !reason.lt.0
! Get the solution and copy into ctl%x_c array
call TaoGetSolutionVector(tao, MyState, ierr)
@@ -183,7 +220,12 @@ end subroutine tao_minimizer
subroutine MyFuncAndGradient(tao, MyState, CostFunc, Grad, dummy, ierr)
+#include
+#if PETSC_VERSION_GE(3,8,0)
+#include "petsc/finclude/petscvec.h"
+#else
#include "petsc/finclude/petscvecdef.h"
+#endif
use set_knd
use drv_str
diff --git a/veof_chl.f90 b/veof_chl.f90
index f6ddbc8..810def3 100644
--- a/veof_chl.f90
+++ b/veof_chl.f90
@@ -33,35 +33,72 @@ subroutine veof_chl
use drv_str
use grd_str
use eof_str
+ use mpi_str
implicit none
- INTEGER(i4) :: i, j, k, l,n, k1
+ INTEGER(i4) :: i, j, k, l,n, k1, my_km, MyNEofs, ierr
REAL(r8), DIMENSION ( grd%im, grd%jm) :: egm
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: eva
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: evc
+ my_km = 0
+ if((drv%chl_assim .eq.1) .and. (drv%multiv .eq. 0)) then
+ MyNEofs = ros%neof_chl
+ my_km = grd%km
+ else if((drv%chl_assim .eq.0) .and. (drv%multiv .eq. 1)) then
+ MyNEofs = ros%neof_multi
+ my_km = ros%kmchl
+ endif
+
+ if(my_km .eq. 0) then
+ if(MyId .eq. 0) then
+ write(drv%dia,*) "Error! my_km for chlorophyll not setted"
+ write(drv%dia,*) "chl_assim e multiv flags should be alternatively valid"
+ write(drv%dia,*) ""
+ write(*,*) "Error! my_km for chlorophyll not setted! Aborting"
+ write(*,*) "chl_assim e multiv flags should be alternatively valid"
+ write(*,*) ""
+ endif
+ call MPI_Barrier(Var3DCommunicator, ierr)
+ call MPI_Abort(Var3DCommunicator,-1,ierr)
+ endif
+
+ ALLOCATE (eva(ros%nreg,MyNEofs)); eva = huge(eva(1,1))
+ ALLOCATE (evc(ros%nreg,my_km,MyNEofs)); evc = huge(evc(1,1,1))
+
+ if((drv%chl_assim .eq.1) .and. (drv%multiv .eq. 0)) then
+ eva(:,:) = ros%eva_chl(:,:)
+ evc(:,:,:) = ros%evc_chl(:,:,:)
+ else if((drv%chl_assim .eq.0) .and. (drv%multiv .eq. 1)) then
+ eva(:,:) = ros%eva_multi(:,:)
+ evc(:,1:my_km,:) = ros%evc_multi(:,1:my_km,:)
+ endif
grd%chl(:,:,:) = 0.0
!cdir noconcur
- do n=1,ros%neof_chl
+ do n=1,MyNEofs
egm(:,:) = 0.0
do j=1,grd%jm
do i=1,grd%im
- egm(i,j) = ros%eva_chl(grd%reg(i,j),n) * grd%ro( i, j, n)
+ egm(i,j) = eva(grd%reg(i,j),n) * grd%ro( i, j, n)
enddo
enddo
! 3D variables
- do k=1,grd%km ! OMP
+ do k=1,my_km ! OMP
k1 = k1 + 1
do j=1,grd%jm
do i=1,grd%im
- grd%chl(i,j,k) = grd%chl(i,j,k) + ros%evc_chl(grd%reg(i,j),k,n) * egm(i,j)
+ grd%chl(i,j,k) = grd%chl(i,j,k) + evc(grd%reg(i,j),k,n) * egm(i,j)
enddo
enddo
enddo
enddo
+
+ DEALLOCATE(eva,evc)
end subroutine veof_chl
diff --git a/veof_chl_ad.f90 b/veof_chl_ad.f90
index aa76325..6f1d8b0 100644
--- a/veof_chl_ad.f90
+++ b/veof_chl_ad.f90
@@ -33,13 +33,49 @@ subroutine veof_chl_ad
use drv_str
use grd_str
use eof_str
+ use mpi_str
implicit none
- INTEGER(i4) :: i, j, k, l, n, k1
+ INTEGER(i4) :: i, j, k, l, n, my_km, MyNEofs, ierr!, k1
REAL(r8), DIMENSION ( grd%im, grd%jm) :: egm
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: eva
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: evc
-do n=1,ros%neof_chl
+my_km = 0
+if((drv%chl_assim .eq.1) .and. (drv%multiv .eq. 0)) then
+ MyNEofs = ros%neof_chl
+ my_km = grd%km
+else if((drv%chl_assim .eq.0) .and. (drv%multiv .eq. 1)) then
+ MyNEofs = ros%neof_multi
+ my_km = ros%kmchl
+endif
+
+if(my_km .eq. 0) then
+ if(MyId .eq. 0) then
+ write(drv%dia,*) "Error! my_km for chlorophyll not setted"
+ write(drv%dia,*) "chl_assim e multiv flags should be alternatively valid"
+ write(drv%dia,*) ""
+ write(*,*) "Error! my_km for chlorophyll not setted! Aborting"
+ write(*,*) "chl_assim e multiv flags should be alternatively valid"
+ write(*,*) ""
+ endif
+ call MPI_Barrier(Var3DCommunicator, ierr)
+ call MPI_Abort(Var3DCommunicator,-1,ierr)
+endif
+
+ALLOCATE (eva(ros%nreg,MyNEofs)); eva = huge(eva(1,1))
+ALLOCATE (evc(ros%nreg,my_km,MyNEofs)); evc = huge(evc(1,1,1))
+
+if((drv%chl_assim .eq.1) .and. (drv%multiv .eq. 0)) then
+ eva(:,:) = ros%eva_chl(:,:)
+ evc(:,:,:) = ros%evc_chl(:,:,:)
+else if((drv%chl_assim .eq.0) .and. (drv%multiv .eq. 1)) then
+ eva(:,:) = ros%eva_multi(:,:)
+ evc(:,1:my_km,:) = ros%evc_multi(:,1:my_km,:)
+endif
+
+do n=1,MyNEofs
grd%ro_ad(:,:,n) = 0.0 ! OMP
enddo
@@ -47,18 +83,18 @@ subroutine veof_chl_ad
!$OMP PRIVATE(i,j,k,k1,n) &
!$OMP PRIVATE(egm)
!$OMP DO
- do n=1,ros%neof_chl
+ do n=1,MyNEofs
egm(:,:) = 0.0
! 3D variables
- k1 = 0
+ ! k1 = 0
- do k=1,grd%km ! OMP
- k1 = k1 + 1
+ do k=1,my_km ! OMP
+ ! k1 = k1 + 1
do j=1,grd%jm
do i=1,grd%im
- egm(i,j) = egm(i,j) + ros%evc_chl(grd%reg(i,j), k,n) * grd%chl_ad(i,j,k)
+ egm(i,j) = egm(i,j) + evc(grd%reg(i,j), k,n) * grd%chl_ad(i,j,k)
enddo
enddo
enddo
@@ -66,7 +102,7 @@ subroutine veof_chl_ad
do j=1,grd%jm
do i=1,grd%im
- egm(i,j) = ros%eva_chl(grd%reg(i,j),n) * egm(i,j)
+ egm(i,j) = eva(grd%reg(i,j),n) * egm(i,j)
enddo
enddo
@@ -85,5 +121,6 @@ subroutine veof_chl_ad
!$OMP END DO
!$OMP END PARALLEL
+DEALLOCATE(eva,evc)
end subroutine veof_chl_ad
diff --git a/veof_multiv_ad.f90 b/veof_multiv_ad.f90
new file mode 100644
index 0000000..7e8f802
--- /dev/null
+++ b/veof_multiv_ad.f90
@@ -0,0 +1,98 @@
+subroutine veof_multiv_ad
+
+!---------------------------------------------------------------------------
+! !
+! Copyright 2006 Srdjan Dobricic, CMCC, Bologna !
+! !
+! This file is part of OceanVar. !
+! !
+! OceanVar is free software: you can redistribute it and/or modify. !
+! it under the terms of the GNU General Public License as published by !
+! the Free Software Foundation, either version 3 of the License, or !
+! (at your option) any later version. !
+! !
+! OceanVar is distributed in the hope that it will be useful, !
+! but WITHOUT ANY WARRANTY; without even the implied warranty of !
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
+! GNU General Public License for more details. !
+! !
+! You should have received a copy of the GNU General Public License !
+! along with OceanVar. If not, see . !
+! !
+!---------------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! !
+! Vertical transformation (adjoint) !
+! !
+! Version 1: S.Dobricic 2006 !
+!-----------------------------------------------------------------------
+
+
+ use set_knd
+ use drv_str
+ use grd_str
+ use eof_str
+ use mpi_str
+
+ implicit none
+
+ INTEGER(i4) :: i, j, k, l, n, my_km !, k1
+ REAL(r8), DIMENSION ( grd%im, grd%jm) :: egm
+
+my_km = ros%kmchl
+
+
+do n=1,ros%neof_multi
+ grd%ro_ad(:,:,n) = 0.0 ! OMP
+enddo
+
+!$OMP PARALLEL &
+!$OMP PRIVATE(i,j,k,k1,n) &
+!$OMP PRIVATE(egm)
+!$OMP DO
+do n=1,ros%neof_multi
+
+ egm(:,:) = 0.0
+
+ ! 3D variables
+ ! k1 = 0
+
+ do k=1,grd%km ! OMP
+ ! k1 = k1 + 1
+ do j=1,grd%jm
+ do i=1,grd%im
+ if(k.le.my_km) then
+ egm(i,j) = egm(i,j) + ros%evc_multi(grd%reg(i,j), k,n) * grd%chl_ad(i,j,k)
+
+ endif
+ egm(i,j) = egm(i,j) + ros%evc_multi(grd%reg(i,j),k+my_km,n) * grd%n3n_ad(i,j,k)
+ enddo
+ enddo
+ enddo
+
+
+
+ do j=1,grd%jm
+ do i=1,grd%im
+ egm(i,j) = ros%eva_multi(grd%reg(i,j),n) * egm(i,j)
+ enddo
+ enddo
+
+ !cdir serial
+ ! 3D variables
+ ! do l=n,ros%neof
+ do j=1,grd%jm
+ do i=1,grd%im
+ grd%ro_ad(i,j,n) = grd%ro_ad(i,j,n) + egm(i,j)
+ enddo
+ enddo
+ ! enddo
+ !cdir end serial
+
+enddo
+!$OMP END DO
+!$OMP END PARALLEL
+
+
+end subroutine veof_multiv_ad
diff --git a/veof_nut.f90 b/veof_nut.f90
index 98b19eb..49fe217 100644
--- a/veof_nut.f90
+++ b/veof_nut.f90
@@ -33,24 +33,77 @@ subroutine veof_nut(NutArray, Var)
use drv_str
use grd_str
use eof_str
+ use mpi_str
implicit none
- INTEGER(i4) :: i, j, k, l,n, k1
+ INTEGER(i4) :: i, j, k, l,n, my_km, ierr
REAL(r8), DIMENSION ( grd%im, grd%jm) :: egm
REAL(r8) :: NutArray(grd%im,grd%jm,grd%km)
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: eva
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: evc
INTEGER(I4) :: MyNEofs, offset
CHARACTER :: Var
NutArray(:,:,:) = 0.0
+ my_km = 0
offset = 0
- if(Var .eq. 'N') then
+ if((drv%nut .eq.1) .and. (drv%multiv .eq. 0)) then
+ my_km = grd%km
+ if(Var .eq. 'N') then
MyNEofs = ros%neof_n3n
offset = ros%neof_chl
- else
+ else
MyNEofs = ros%neof_o2o
offset = ros%neof_chl + ros%neof_n3n
+ endif
+ else if((drv%nut .eq.0) .and. (drv%multiv .eq. 1)) then
+ if(Var .eq. 'N') then
+ my_km = grd%km
+ MyNEofs = ros%neof_multi
+ else
+ if(MyId .eq. 0) then
+ write(drv%dia,*) "Error! Only nitrate multivariate assimilation implemented"
+ write(drv%dia,*) ""
+ write(*,*) "Error! Only nitrate multivariate assimilation implemented! Aborting"
+ write(*,*) ""
+ endif
+ call MPI_Barrier(Var3DCommunicator, ierr)
+ call MPI_Abort(Var3DCommunicator,-1,ierr)
+ endif
+ endif
+
+
+ if(my_km .eq. 0) then
+ if(MyId .eq. 0) then
+ write(drv%dia,*) "Error! my_km for nutrient not setted"
+ write(drv%dia,*) "drv%nut e multiv flags should be alternatively valid"
+ write(drv%dia,*) ""
+ write(*,*) "Error! my_km for nutrient not setted! Aborting"
+ write(*,*) "drv%nut e multiv flags should be alternatively valid"
+ write(*,*) ""
+ endif
+ call MPI_Barrier(Var3DCommunicator, ierr)
+ call MPI_Abort(Var3DCommunicator,-1,ierr)
+ endif
+
+ ALLOCATE (eva(ros%nreg,MyNEofs)); eva = huge(eva(1,1))
+ ALLOCATE (evc(ros%nreg,my_km,MyNEofs)); evc = huge(evc(1,1,1))
+
+ if((drv%nut .eq.1) .and. (drv%multiv .eq. 0)) then
+ if(Var .eq. 'N') then
+ eva = ros%eva_n3n
+ evc = ros%evc_n3n
+ else
+ eva = ros%eva_o2o
+ evc = ros%evc_o2o
+ endif
+ else if((drv%nut .eq.0) .and. (drv%multiv .eq. 1)) then
+ if(Var .eq. 'N') then
+ eva = ros%eva_multi
+ evc(:,1:my_km,:) = ros%evc_multi(:,ros%kmchl+1:ros%kmchl+grd%km,:)
+ endif
endif
!cdir noconcur
@@ -60,27 +113,19 @@ subroutine veof_nut(NutArray, Var)
do j=1,grd%jm
do i=1,grd%im
- if(Var .eq. 'N') then
- egm(i,j) = ros%eva_n3n(grd%reg(i,j),n) * grd%ro( i, j, n+offset)
- else
- egm(i,j) = ros%eva_o2o(grd%reg(i,j),n) * grd%ro( i, j, n+offset)
- endif
+ egm(i,j) = eva(grd%reg(i,j),n) * grd%ro( i, j, n+offset)
enddo
enddo
! 3D variables
- do k=1,grd%km ! OMP
- k1 = k1 + 1
+ do k=1,my_km ! OMP
do j=1,grd%jm
do i=1,grd%im
- if(Var .eq. 'N') then
- NutArray(i,j,k) = NutArray(i,j,k) + ros%evc_n3n(grd%reg(i,j),k,n) * egm(i,j)
- else
- NutArray(i,j,k) = NutArray(i,j,k) + ros%evc_o2o(grd%reg(i,j),k,n) * egm(i,j)
- endif
+ NutArray(i,j,k) = NutArray(i,j,k) + evc(grd%reg(i,j),k,n) * egm(i,j)
enddo
enddo
enddo
enddo
-
+
+DEALLOCATE(eva,evc)
end subroutine veof_nut
diff --git a/veof_nut_ad.f90 b/veof_nut_ad.f90
index cc55ca7..858a8a5 100644
--- a/veof_nut_ad.f90
+++ b/veof_nut_ad.f90
@@ -36,19 +36,51 @@ subroutine veof_nut_ad(NutArrayAd, Var)
implicit none
- INTEGER(i4) :: i, j, k, l, n, k1, offset
+ INTEGER(i4) :: i, j, k, l, n, offset, my_km
REAL(r8), DIMENSION ( grd%im, grd%jm) :: egm
REAL(r8) :: NutArrayAd(grd%im,grd%jm,grd%km)
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: eva
+ REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: evc
CHARACTER :: Var
INTEGER :: MyNEofs
+ my_km = 0
+ ! Altrove usato grd%km come limite per assimilazione nit qui ro%kmnit
+ ! Da correggere o fare un check
offset = 0
- if(Var .eq. 'N') then
- MyNEofs = ros%neof_n3n
- offset = ros%neof_chl
- else
- MyNEofs = ros%neof_o2o
- offset = ros%neof_chl + ros%neof_n3n
+ if((drv%nut .eq.1) .and. (drv%multiv .eq. 0)) then
+ my_km = grd%km
+ if(Var .eq. 'N') then
+ MyNEofs = ros%neof_n3n
+ offset = ros%neof_chl
+ else
+ MyNEofs = ros%neof_o2o
+ offset = ros%neof_chl + ros%neof_n3n
+ endif
+ else if((drv%nut .eq.0) .and. (drv%multiv .eq. 1)) then
+ if(Var .eq. 'N') then
+ my_km = ros%kmnit
+ MyNEofs = ros%neof_multi
+ endif
+ endif
+
+
+ ALLOCATE (eva(ros%nreg,MyNEofs)); eva = huge(eva(1,1))
+ ALLOCATE (evc(ros%nreg,my_km,MyNEofs)); evc = huge(evc(1,1,1))
+
+ if((drv%nut .eq.1) .and. (drv%multiv .eq. 0)) then
+ if(Var .eq. 'N') then
+ eva = ros%eva_n3n
+ evc = ros%evc_n3n
+ else
+ eva = ros%eva_o2o
+ evc = ros%evc_o2o
+ endif
+ else if((drv%nut .eq.0) .and. (drv%multiv .eq. 1)) then
+ if(Var .eq. 'N') then
+ eva = ros%eva_multi
+ evc(:,1:my_km,:) = ros%evc_multi(:,ros%kmchl+1:ros%kmchl+ros%kmnit,:)
+ endif
endif
do n=1,MyNEofs
@@ -64,17 +96,11 @@ subroutine veof_nut_ad(NutArrayAd, Var)
egm(:,:) = 0.0
! 3D variables
- k1 = 0
- do k=1,grd%km ! OMP
- k1 = k1 + 1
+ do k=1,my_km ! OMP
do j=1,grd%jm
do i=1,grd%im
- if(Var .eq. 'N') then
- egm(i,j) = egm(i,j) + ros%evc_n3n(grd%reg(i,j), k,n) * NutArrayAd(i,j,k)
- else
- egm(i,j) = egm(i,j) + ros%evc_o2o(grd%reg(i,j), k,n) * NutArrayAd(i,j,k)
- endif
+ egm(i,j) = egm(i,j) + evc(grd%reg(i,j), k,n) * NutArrayAd(i,j,k)
enddo
enddo
enddo
@@ -82,14 +108,10 @@ subroutine veof_nut_ad(NutArrayAd, Var)
do j=1,grd%jm
do i=1,grd%im
- if(Var .eq. 'N') then
- egm(i,j) = ros%eva_n3n(grd%reg(i,j),n) * egm(i,j)
- else
- egm(i,j) = ros%eva_o2o(grd%reg(i,j),n) * egm(i,j)
- endif
+ egm(i,j) = eva(grd%reg(i,j),n) * egm(i,j)
enddo
enddo
-
+
do j=1,grd%jm
do i=1,grd%im
grd%ro_ad(i,j,n+offset) = grd%ro_ad(i,j,n+offset) + egm(i,j)
@@ -100,5 +122,6 @@ subroutine veof_nut_ad(NutArrayAd, Var)
!$OMP END DO
!$OMP END PARALLEL
+DEALLOCATE(eva,evc)
end subroutine veof_nut_ad
diff --git a/ver_hor_chl.f90 b/ver_hor_chl.f90
index 539d18a..6ed94b5 100644
--- a/ver_hor_chl.f90
+++ b/ver_hor_chl.f90
@@ -44,7 +44,7 @@ subroutine ver_hor_chl
implicit none
- INTEGER(i4) :: i,j,k, ione, l
+ INTEGER(i4) :: i,j,k, ione, l, my_km, myimax, myjmax
INTEGER :: jp, SurfaceIndex, TmpOffset, LinearIndex
INTEGER(i4) :: iProc, ierr
type(DoubleGrid), allocatable, dimension(:,:,:) :: SendBuf3D
@@ -59,9 +59,17 @@ subroutine ver_hor_chl
!return
! goto 103 !No Vh
+ my_km = grd%km
+ myimax = grd%imax
+ myjmax = grd%jmax
+ if(drv%multiv.eq.1) then
+ my_km = ros%kmchl
+ myimax = grd%imax_chl
+ myjmax = grd%jmax_chl
+ endif
! ---
! Load temporary arrays
- do k=1,grd%km
+ do k=1,my_km
grd%chl_ad(:,:,k) = grd%chl(:,:,k)
enddo
@@ -73,46 +81,46 @@ subroutine ver_hor_chl
! y direction
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
grd%chl_ad(:,:,k) = grd%chl_ad(:,:,k) * grd%scy(:,:,k)
enddo
! Apply recursive filter in y direction
- call rcfl_y_ad( localRow, GlobalCol, grd%km, grd%jmax, grd%aey, grd%bey, grd%chl_ad, grd%jnx, grd%jmx)
-
+ call rcfl_y_ad( localRow, GlobalCol, my_km, myjmax, grd%aey(:,:,1:my_km), &
+ grd%bey(:,:,1:my_km), grd%chl_ad, grd%jnx(:,:,1:my_km), grd%jmx(1:my_km))
! ---
! x direction
if(NumProcI .gt. 1) then
- ALLOCATE(SendBuf3D(grd%km, grd%im, grd%jm))
- ALLOCATE( RecBuf1D(grd%km*GlobalRow*localCol))
- ALLOCATE( DefBufChl(GlobalRow, localCol, grd%km))
- ALLOCATE( DefBufChlAd(GlobalRow, localCol, grd%km))
+ ALLOCATE(SendBuf3D(my_km, grd%im, grd%jm))
+ ALLOCATE( RecBuf1D(my_km*GlobalRow*localCol))
+ ALLOCATE( DefBufChl(GlobalRow, localCol, my_km))
+ ALLOCATE( DefBufChlAd(GlobalRow, localCol, my_km))
- do k=1,grd%km
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
SendBuf3D(k,i,j)%chl = grd%chl(i,j,k)
end do
end do
end do
- do k=1,grd%km
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
SendBuf3D(k,i,j)%chl_ad = grd%chl_ad(i,j,k)
end do
end do
end do
+
+ call MPI_Alltoallv(SendBuf3D, SendCountX3D_chl, SendDisplX3D_chl, MyPair, &
+ RecBuf1D, RecCountX3D_chl, RecDisplX3D_chl, MyPair, Var3DCommunicator, ierr)
- call MPI_Alltoallv(SendBuf3D, SendCountX3D, SendDisplX3D, MyPair, &
- RecBuf1D, RecCountX3D, RecDisplX3D, MyPair, Var3DCommunicator, ierr)
-
- SurfaceIndex = localCol*grd%km
+ SurfaceIndex = localCol*my_km
do j=1,localCol
do iProc=0, NumProcI-1
- TmpOffset = RecDisplX3D(iProc+1)/SurfaceIndex
- do i=1,RecCountX3D(iProc+1)/SurfaceIndex
- LinearIndex = (i-1)*grd%km + (j-1)*RecCountX3D(iProc+1)/localCol + RecDisplX3D(iProc+1)
- do k=1,grd%km
+ TmpOffset = RecDisplX3D_chl(iProc+1)/SurfaceIndex
+ do i=1,RecCountX3D_chl(iProc+1)/SurfaceIndex
+ LinearIndex = (i-1)*my_km + (j-1)*RecCountX3D_chl(iProc+1)/localCol + RecDisplX3D_chl(iProc+1)
+ do k=1,my_km
DefBufChl(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl
end do
end do
@@ -121,10 +129,10 @@ subroutine ver_hor_chl
end do
do j=1,localCol
do iProc=0, NumProcI-1
- TmpOffset = RecDisplX3D(iProc+1)/SurfaceIndex
- do i=1,RecCountX3D(iProc+1)/SurfaceIndex
- LinearIndex = (i-1)*grd%km + (j-1)*RecCountX3D(iProc+1)/localCol + RecDisplX3D(iProc+1)
- do k=1,grd%km
+ TmpOffset = RecDisplX3D_chl(iProc+1)/SurfaceIndex
+ do i=1,RecCountX3D_chl(iProc+1)/SurfaceIndex
+ LinearIndex = (i-1)*my_km + (j-1)*RecCountX3D_chl(iProc+1)/localCol + RecDisplX3D_chl(iProc+1)
+ do k=1,my_km
DefBufChlAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl_ad
end do
end do
@@ -133,20 +141,22 @@ subroutine ver_hor_chl
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
DefBufChlAd(:,:,k) = DefBufChlAd(:,:,k) * grd%scx(:,:,k)
enddo
- call rcfl_x_ad( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, DefBufChlAd, grd%inx, grd%imx)
+ call rcfl_x_ad( GlobalRow, localCol, my_km, myimax, grd%aex(:,:,1:my_km), &
+ grd%bex(:,:,1:my_km), DefBufChlAd, grd%inx(:,:,1:my_km), grd%imx(1:my_km))
else
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
grd%chl_ad(:,:,k) = grd%chl_ad(:,:,k) * grd%scx(:,:,k)
enddo
- call rcfl_x_ad( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, grd%chl_ad, grd%inx, grd%imx)
+ call rcfl_x_ad( GlobalRow, localCol, my_km, myimax, grd%aex(:,:,1:my_km), &
+ grd%bex(:,:,1:my_km), grd%chl_ad, grd%inx(:,:,1:my_km), grd%imx(1:my_km))
end if
@@ -156,25 +166,26 @@ subroutine ver_hor_chl
! x direction
if(NumProcI .gt. 1) then
- call rcfl_x( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, DefBufChl, grd%inx, grd%imx)
+ call rcfl_x( GlobalRow, localCol, my_km, myimax, grd%aex(:,:,1:my_km), &
+ grd%bex(:,:,1:my_km), DefBufChl, grd%inx(:,:,1:my_km), grd%imx(1:my_km))
- do k=1,grd%km
+ do k=1,my_km
DefBufChl(:,:,k) = DefBufChl(:,:,k) * grd%scx(:,:,k)
enddo
! Reordering data to send back
DEALLOCATE(SendBuf3D, RecBuf1D)
- ALLOCATE(SendBuf3D(grd%km, localCol, GlobalRow))
- ALLOCATE( RecBuf1D(grd%km*grd%jm*grd%im))
+ ALLOCATE(SendBuf3D(my_km, localCol, GlobalRow))
+ ALLOCATE( RecBuf1D(my_km*grd%jm*grd%im))
- do k=1,grd%km
+ do k=1,my_km
do j=1,localCol
do i=1,GlobalRow
SendBuf3D(k,j,i)%chl = DefBufChl(i,j,k)
end do
end do
end do
- do k=1,grd%km
+ do k=1,my_km
do j=1,localCol
do i=1,GlobalRow
SendBuf3D(k,j,i)%chl_ad = DefBufChlAd(i,j,k)
@@ -182,16 +193,16 @@ subroutine ver_hor_chl
end do
end do
- call MPI_Alltoallv(SendBuf3D, RecCountX3D, RecDisplX3D, MyPair, &
- RecBuf1D, SendCountX3D, SendDisplX3D, MyPair, Var3DCommunicator, ierr)
+ call MPI_Alltoallv(SendBuf3D, RecCountX3D_chl, RecDisplX3D_chl, MyPair, &
+ RecBuf1D, SendCountX3D_chl, SendDisplX3D_chl, MyPair, Var3DCommunicator, ierr)
- SurfaceIndex = grd%im*grd%km
+ SurfaceIndex = grd%im*my_km
do i=1,grd%im
do iProc=0, NumProcI-1
- TmpOffset = SendDisplX3D(iProc+1)/SurfaceIndex
- do j=1,SendCountX3D(iProc+1)/SurfaceIndex
- LinearIndex = (j-1)*grd%km +(i-1)*SendCountX3D(iProc+1)/grd%im + SendDisplX3D(iProc+1)
- do k=1,grd%km
+ TmpOffset = SendDisplX3D_chl(iProc+1)/SurfaceIndex
+ do j=1,SendCountX3D_chl(iProc+1)/SurfaceIndex
+ LinearIndex = (j-1)*my_km +(i-1)*SendCountX3D_chl(iProc+1)/grd%im + SendDisplX3D_chl(iProc+1)
+ do k=1,my_km
grd%chl(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl
end do
end do
@@ -199,10 +210,10 @@ subroutine ver_hor_chl
end do
do i=1,grd%im
do iProc=0, NumProcI-1
- TmpOffset = SendDisplX3D(iProc+1)/SurfaceIndex
- do j=1,SendCountX3D(iProc+1)/SurfaceIndex
- LinearIndex = (j-1)*grd%km +(i-1)*SendCountX3D(iProc+1)/grd%im + SendDisplX3D(iProc+1)
- do k=1,grd%km
+ TmpOffset = SendDisplX3D_chl(iProc+1)/SurfaceIndex
+ do j=1,SendCountX3D_chl(iProc+1)/SurfaceIndex
+ LinearIndex = (j-1)*my_km +(i-1)*SendCountX3D_chl(iProc+1)/grd%im + SendDisplX3D_chl(iProc+1)
+ do k=1,my_km
grd%chl_ad(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl_ad
end do
end do
@@ -213,9 +224,10 @@ subroutine ver_hor_chl
else ! NumProcI .eq. 1
- call rcfl_x( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, grd%chl, grd%inx, grd%imx)
+ call rcfl_x( GlobalRow, localCol, my_km, myimax, grd%aex(:,:,1:my_km), &
+ grd%bex(:,:,1:my_km), grd%chl, grd%inx(:,:,1:my_km), grd%imx(1:my_km))
- do k=1,grd%km
+ do k=1,my_km
grd%chl(:,:,k) = grd%chl(:,:,k) * grd%scx(:,:,k)
enddo
@@ -224,23 +236,24 @@ subroutine ver_hor_chl
! ---
! y direction
! Apply recursive filter in y direction
- call rcfl_y( localRow, GlobalCol, grd%km, grd%jmax, grd%aey, grd%bey, grd%chl, grd%jnx, grd%jmx)
+ call rcfl_y( localRow, GlobalCol, my_km, myjmax, grd%aey(:,:,1:my_km), &
+ grd%bey(:,:,1:my_km), grd%chl, grd%jnx(:,:,1:my_km), grd%jmx(1:my_km))
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
grd%chl(:,:,k) = grd%chl(:,:,k) * grd%scy(:,:,k)
enddo
! ---
! Average
- do k=1,grd%km
+ do k=1,my_km
grd%chl(:,:,k) = (grd%chl(:,:,k) + grd%chl_ad(:,:,k) ) * 0.5
enddo
! ---
! Scale for boundaries
- do k=1,grd%km
+ do k=1,my_km
grd%chl(:,:,k) = grd%chl(:,:,k) * grd%msk(:,:,k)
enddo
diff --git a/ver_hor_chl_ad.f90 b/ver_hor_chl_ad.f90
index 3eb3a53..6205412 100644
--- a/ver_hor_chl_ad.f90
+++ b/ver_hor_chl_ad.f90
@@ -23,7 +23,7 @@ subroutine ver_hor_chl_ad
implicit none
- INTEGER(i4) :: i,j,k, ione, l
+ INTEGER(i4) :: i,j,k, ione, l, my_km, myimax, myjmax
INTEGER :: jp, SurfaceIndex, TmpOffset, LinearIndex
INTEGER(i4) :: iProc, ierr
type(DoubleGrid), allocatable, dimension(:,:,:) :: SendBuf3D
@@ -34,16 +34,24 @@ subroutine ver_hor_chl_ad
! goto 103 ! No Vh
+ my_km = grd%km
+ myimax = grd%imax
+ myjmax = grd%jmax
+ if(drv%multiv.eq.1) then
+ my_km = ros%kmchl
+ myimax = grd%imax_chl
+ myjmax = grd%jmax_chl
+ endif
! ---
! Scale for boundaries
- do k=1,grd%km
+ do k=1,my_km
grd%chl_ad(:,:,k) = grd%chl_ad(:,:,k) * grd%msk(:,:,k)
enddo
! ---
! Load temporary arrays
- do k=1,grd%km
+ do k=1,my_km
grd%chl(:,:,k) = grd%chl_ad(:,:,k)
enddo
@@ -51,29 +59,30 @@ subroutine ver_hor_chl_ad
! y direction
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
grd%chl_ad(:,:,k) = grd%chl_ad(:,:,k) * grd%scy(:,:,k)
enddo
! Apply recursive filter in y direction
- call rcfl_y_ad( localRow, GlobalCol, grd%km, grd%jmax, grd%aey, grd%bey, grd%chl_ad, grd%jnx, grd%jmx)
+ call rcfl_y_ad( localRow, GlobalCol, my_km, myjmax, grd%aey(:,:,1:my_km), &
+ grd%bey(:,:,1:my_km), grd%chl_ad, grd%jnx(:,:,1:my_km), grd%jmx(1:my_km))
! ---
! x direction
if(NumProcI .gt. 1) then
- ALLOCATE(SendBuf3D(grd%km, grd%im, grd%jm))
- ALLOCATE( RecBuf1D(grd%km*localCol*GlobalRow))
- ALLOCATE(DefBufChl(GlobalRow, localCol, grd%km))
- ALLOCATE(DefBufChlAd(GlobalRow, localCol, grd%km))
+ ALLOCATE(SendBuf3D(my_km, grd%im, grd%jm))
+ ALLOCATE( RecBuf1D(my_km*localCol*GlobalRow))
+ ALLOCATE(DefBufChl(GlobalRow, localCol, my_km))
+ ALLOCATE(DefBufChlAd(GlobalRow, localCol, my_km))
- do k=1,grd%km
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
SendBuf3D(k,i,j)%chl = grd%chl(i,j,k)
end do
end do
end do
- do k=1,grd%km
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
SendBuf3D(k,i,j)%chl_ad = grd%chl_ad(i,j,k)
@@ -81,16 +90,16 @@ subroutine ver_hor_chl_ad
end do
end do
- call MPI_Alltoallv(SendBuf3D, SendCountX3D, SendDisplX3D, MyPair, &
- RecBuf1D, RecCountX3D, RecDisplX3D, MyPair, Var3DCommunicator, ierr)
+ call MPI_Alltoallv(SendBuf3D, SendCountX3D_chl, SendDisplX3D_chl, MyPair, &
+ RecBuf1D, RecCountX3D_chl, RecDisplX3D_chl, MyPair, Var3DCommunicator, ierr)
- SurfaceIndex = localCol*grd%km
+ SurfaceIndex = localCol*my_km
do j=1,localCol
do iProc=0, NumProcI-1
- TmpOffset = RecDisplX3D(iProc+1)/SurfaceIndex
- do i=1,RecCountX3D(iProc+1)/SurfaceIndex
- LinearIndex = (i-1)*grd%km + (j-1)*RecCountX3D(iProc+1)/localCol + RecDisplX3D(iProc+1)
- do k=1,grd%km
+ TmpOffset = RecDisplX3D_chl(iProc+1)/SurfaceIndex
+ do i=1,RecCountX3D_chl(iProc+1)/SurfaceIndex
+ LinearIndex = (i-1)*my_km + (j-1)*RecCountX3D_chl(iProc+1)/localCol + RecDisplX3D_chl(iProc+1)
+ do k=1,my_km
DefBufChl(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl
end do
end do
@@ -98,10 +107,10 @@ subroutine ver_hor_chl_ad
end do
do j=1,localCol
do iProc=0, NumProcI-1
- TmpOffset = RecDisplX3D(iProc+1)/SurfaceIndex
- do i=1,RecCountX3D(iProc+1)/SurfaceIndex
- LinearIndex = (i-1)*grd%km + (j-1)*RecCountX3D(iProc+1)/localCol + RecDisplX3D(iProc+1)
- do k=1,grd%km
+ TmpOffset = RecDisplX3D_chl(iProc+1)/SurfaceIndex
+ do i=1,RecCountX3D_chl(iProc+1)/SurfaceIndex
+ LinearIndex = (i-1)*my_km + (j-1)*RecCountX3D_chl(iProc+1)/localCol + RecDisplX3D_chl(iProc+1)
+ do k=1,my_km
DefBufChlAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl_ad
end do
end do
@@ -110,20 +119,22 @@ subroutine ver_hor_chl_ad
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
DefBufChlAd(:,:,k) = DefBufChlAd(:,:,k) * grd%scx(:,:,k)
enddo
- call rcfl_x_ad( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, DefBufChlAd, grd%inx, grd%imx)
+ call rcfl_x_ad( GlobalRow, localCol, my_km, myimax, grd%aex(:,:,1:my_km), &
+ grd%bex(:,:,1:my_km), DefBufChlAd, grd%inx(:,:,1:my_km), grd%imx(1:my_km))
else ! NumProcI .eq. 1
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
grd%chl_ad(:,:,k) = grd%chl_ad(:,:,k) * grd%scx(:,:,k)
enddo
- call rcfl_x_ad( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, grd%chl_ad, grd%inx, grd%imx)
+ call rcfl_x_ad( GlobalRow, localCol, my_km, myimax, grd%aex(:,:,1:my_km), &
+ grd%bex(:,:,1:my_km), grd%chl_ad, grd%inx(:,:,1:my_km), grd%imx(1:my_km))
end if
@@ -131,27 +142,28 @@ subroutine ver_hor_chl_ad
! x direction
if(NumProcI .gt. 1) then
- call rcfl_x( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, DefBufChl, grd%inx, grd%imx)
+ call rcfl_x( GlobalRow, localCol, my_km, myimax, grd%aex(:,:,1:my_km), &
+ grd%bex(:,:,1:my_km), DefBufChl, grd%inx(:,:,1:my_km), grd%imx(1:my_km))
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
DefBufChl(:,:,k) = DefBufChl(:,:,k) * grd%scx(:,:,k)
enddo
! Reordering data to send back
DEALLOCATE(SendBuf3D, RecBuf1D)
- ALLOCATE(SendBuf3D(grd%km, localCol, GlobalRow))
- ALLOCATE( RecBuf1D(grd%km*grd%jm*grd%im))
+ ALLOCATE(SendBuf3D(my_km, localCol, GlobalRow))
+ ALLOCATE( RecBuf1D(my_km*grd%jm*grd%im))
- do k=1,grd%km
+ do k=1,my_km
do j=1,localCol
do i=1,GlobalRow
SendBuf3D(k,j,i)%chl = DefBufChl(i,j,k)
end do
end do
end do
- do k=1,grd%km
+ do k=1,my_km
do j=1,localCol
do i=1,GlobalRow
SendBuf3D(k,j,i)%chl_ad = DefBufChlAd(i,j,k)
@@ -159,16 +171,16 @@ subroutine ver_hor_chl_ad
end do
end do
- call MPI_Alltoallv(SendBuf3D, RecCountX3D, RecDisplX3D, MyPair, &
- RecBuf1D, SendCountX3D, SendDisplX3D, MyPair, Var3DCommunicator, ierr)
+ call MPI_Alltoallv(SendBuf3D, RecCountX3D_chl, RecDisplX3D_chl, MyPair, &
+ RecBuf1D, SendCountX3D_chl, SendDisplX3D_chl, MyPair, Var3DCommunicator, ierr)
- SurfaceIndex = grd%im*grd%km
+ SurfaceIndex = grd%im*my_km
do i=1,grd%im
do iProc=0, NumProcI-1
- TmpOffset = SendDisplX3D(iProc+1)/SurfaceIndex
- do j=1,SendCountX3D(iProc+1)/SurfaceIndex
- LinearIndex = (j-1)*grd%km +(i-1)*SendCountX3D(iProc+1)/grd%im + SendDisplX3D(iProc+1)
- do k=1,grd%km
+ TmpOffset = SendDisplX3D_chl(iProc+1)/SurfaceIndex
+ do j=1,SendCountX3D_chl(iProc+1)/SurfaceIndex
+ LinearIndex = (j-1)*my_km +(i-1)*SendCountX3D_chl(iProc+1)/grd%im + SendDisplX3D_chl(iProc+1)
+ do k=1,my_km
grd%chl(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl
end do
end do
@@ -176,10 +188,10 @@ subroutine ver_hor_chl_ad
end do
do i=1,grd%im
do iProc=0, NumProcI-1
- TmpOffset = SendDisplX3D(iProc+1)/SurfaceIndex
- do j=1,SendCountX3D(iProc+1)/SurfaceIndex
- LinearIndex = (j-1)*grd%km +(i-1)*SendCountX3D(iProc+1)/grd%im + SendDisplX3D(iProc+1)
- do k=1,grd%km
+ TmpOffset = SendDisplX3D_chl(iProc+1)/SurfaceIndex
+ do j=1,SendCountX3D_chl(iProc+1)/SurfaceIndex
+ LinearIndex = (j-1)*my_km +(i-1)*SendCountX3D_chl(iProc+1)/grd%im + SendDisplX3D_chl(iProc+1)
+ do k=1,my_km
grd%chl_ad(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl_ad
end do
end do
@@ -189,11 +201,12 @@ subroutine ver_hor_chl_ad
DEALLOCATE(SendBuf3D, RecBuf1D, DefBufChl, DefBufChlAd)
else ! NumProcI .eq. 1
- call rcfl_x( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, grd%chl, grd%inx, grd%imx)
+ call rcfl_x( GlobalRow, localCol, my_km, myimax, grd%aex(:,:,1:my_km), &
+ grd%bex(:,:,1:my_km), grd%chl, grd%inx(:,:,1:my_km), grd%imx(1:my_km))
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
grd%chl(:,:,k) = grd%chl(:,:,k) * grd%scx(:,:,k)
enddo
end if
@@ -202,18 +215,19 @@ subroutine ver_hor_chl_ad
! ! ---
! ! y direction
! Apply recursive filter in y direction
- call rcfl_y( localRow, GlobalCol, grd%km, grd%jmax, grd%aey, grd%bey, grd%chl, grd%jnx, grd%jmx)
+ call rcfl_y( localRow, GlobalCol, my_km, myjmax, grd%aey(:,:,1:my_km), &
+ grd%bey(:,:,1:my_km), grd%chl, grd%jnx(:,:,1:my_km), grd%jmx(1:my_km))
! ---
! Scale by the scaling factor
- do k=1,grd%km
+ do k=1,my_km
grd%chl(:,:,k) = grd%chl(:,:,k) * grd%scy(:,:,k)
enddo
! ---
! Average
- do k=1,grd%km
+ do k=1,my_km
grd%chl_ad(:,:,k) = (grd%chl_ad(:,:,k) + grd%chl(:,:,k) ) * 0.5
enddo
@@ -221,6 +235,7 @@ subroutine ver_hor_chl_ad
! 103 continue
! ---
! Vertical EOFs
- call veof_chl_ad
+ if(drv%multiv.eq.0) &
+ call veof_chl_ad
end subroutine ver_hor_chl_ad
diff --git a/ver_hor_nut_ad.f90 b/ver_hor_nut_ad.f90
index 93d1872..8e02161 100644
--- a/ver_hor_nut_ad.f90
+++ b/ver_hor_nut_ad.f90
@@ -223,6 +223,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var)
! 103 continue
! ---
! Vertical EOFs
- call veof_nut_ad(NutArrayAd, Var)
+ if(drv%multiv.eq.0) &
+ call veof_nut_ad(NutArrayAd, Var)
end subroutine ver_hor_nut_ad
diff --git a/wrt_chl_stat.f90 b/wrt_chl_stat.f90
index 334748c..4c21db6 100644
--- a/wrt_chl_stat.f90
+++ b/wrt_chl_stat.f90
@@ -2,6 +2,7 @@ subroutine wrt_chl_stat
use set_knd
use grd_str
+ use eof_str
use drv_str
use mpi_str
use bio_str
@@ -10,20 +11,29 @@ subroutine wrt_chl_stat
implicit none
- INTEGER(i4) :: ncid, ierr, i, j, k, l, m, mm
- INTEGER(i4) :: idP, iVar
+ INTEGER(i4) :: ncid, ierr, i, j, k, l, m, mm, my_km
+ INTEGER(i4) :: idP, iVar, idL
INTEGER(I4) :: xid,yid,depid,timeId, idTim
INTEGER :: system, SysErr
INTEGER(kind=MPI_OFFSET_KIND) :: global_im, global_jm, global_km, MyTime
INTEGER(KIND=MPI_OFFSET_KIND) :: MyCountSingle(1), MyStartSingle(1)
- CHARACTER(LEN=37) :: BioRestart
- CHARACTER(LEN=39) :: BioRestartLong
+ CHARACTER(LEN=45) :: BioRestart
+ CHARACTER(LEN=47) :: BioRestartLong
CHARACTER(LEN=6) :: MyVarName
+ CHARACTER(LEN=16) :: LimVarName
+ CHARACTER(LEN=49) :: LimCorrfile
LOGICAL, ALLOCATABLE :: MyConditions(:,:,:,:)
+ INTEGER, ALLOCATABLE :: LimitCorr(:,:,:)
+
+
+ real(r8) :: TmpVal, MyCorr, MyRatio, SMALL
+ real(r4), allocatable, dimension(:,:,:) :: ValuesToTest
+
+ ! bug fix Intel 2018
+ real(r4), allocatable, dimension(:,:,:,:) :: DumpBio
+ integer(KIND=MPI_OFFSET_KIND) :: MyStart_4d(4), MyCount_4d(4)
- real(r8) :: TmpVal, MyCorr, MyRatio,SMALL
- real(r4), allocatable, dimension(:,:,:) :: DumpBio, ValuesToTest
real(r8) :: TimeArr(1)
real(r4) :: MAX_N_CHL, MAX_P_CHL, MAX_P_C, MAX_N_C
real(r4) :: OPT_N_C, OPT_P_C, OPT_S_C, LIM_THETA
@@ -35,12 +45,22 @@ subroutine wrt_chl_stat
MAX_N_C = 1.26e-2*2 ! values from BFMconsortium parametrs document (P.Lazzari)
OPT_N_C = 1.26e-2
OPT_S_C = 0.01 ! values from BFMconsortium parametrs document (P.Lazzari)
- LIM_THETA = 0.01
+ LIM_THETA = 0.001
SMALL = 1.e-5
-
- ALLOCATE(DumpBio(grd%im,grd%jm,grd%km)); DumpBio(:,:,:) = 1.e20
- ALLOCATE(ValuesToTest(grd%im,grd%jm,grd%km)); ValuesToTest(:,:,:) = dble(0.)
- ALLOCATE(MyConditions(grd%im,grd%jm,grd%km,bio%nphy))
+
+ MyStart_4d(1:3) = MyStart(:)
+ MyStart_4d(4) = 1
+ MyCount_4d(1:3) = MyCount(:)
+ MyCount_4d(4) = 1
+
+ my_km = grd%km
+ if(drv%multiv.eq.1) &
+ my_km = ros%kmchl
+
+ ALLOCATE(DumpBio(grd%im,grd%jm,grd%km,1)); DumpBio(:,:,:,:) = 1.e20
+ ALLOCATE(ValuesToTest(grd%im,grd%jm,my_km)); ValuesToTest(:,:,:) = dble(0.)
+ ALLOCATE(MyConditions(grd%im,grd%jm,my_km,bio%nphy))
+ ALLOCATE(LimitCorr(grd%im,grd%jm,grd%km)); LimitCorr(:,:,:) = -99
if(MyId .eq. 0) then
write(drv%dia,*) 'writing chl structure'
@@ -56,7 +76,9 @@ subroutine wrt_chl_stat
MyStartSingle(1) = 1
TimeArr(1) = DA_JulianDate
- do k=1,grd%km
+ LimitCorr(:,:,1:my_km) = 0
+
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
@@ -67,7 +89,11 @@ subroutine wrt_chl_stat
ValuesToTest(i,j,k) = bio%InitialChl(i,j,k) + grd%chl(i,j,k)
if(bio%ApplyConditions) then
if(ValuesToTest(i,j,k) .gt. 10*bio%InitialChl(i,j,k)) then
-
+ if(ValuesToTest(i,j,1) .gt. 0) then
+ LimitCorr(i,j,k) = 1
+ ! print*, ValuesToTest(i,j,1), bio%InitialChl(i,j,1), grd%chl(i,j,1),i,j
+ endif
+
do m=1,bio%ncmp
do l=1,bio%nphy
bio%phy(i,j,k,l,m) = 9.*bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
@@ -97,8 +123,8 @@ subroutine wrt_chl_stat
if(iVar .gt. NPhytoVar) CYCLE
- BioRestart = 'RESTARTS/RST.'//ShortDate//'.'//DA_VarList(iVar)//'.nc'
- BioRestartLong = 'RESTARTS/RST.'//DA_DATE//'.'//DA_VarList(iVar)//'.nc'
+ BioRestart = 'DA__FREQ_1/RST_after.'//ShortDate//'.'//DA_VarList(iVar)//'.nc'
+ BioRestartLong = 'DA__FREQ_1/RST_after.'//DA_DATE//'.'//DA_VarList(iVar)//'.nc'
if(drv%Verbose .eq. 1 .and. MyId .eq. 0) &
print*, "Writing Phyto Restart ", BioRestart
@@ -128,7 +154,7 @@ subroutine wrt_chl_stat
ierr = nf90mpi_enddef(ncid)
if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_enddef'//DA_VarList(iVar), ierr)
- do k=1,grd%km
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
@@ -145,7 +171,8 @@ subroutine wrt_chl_stat
if(TmpVal.gt.SMALL) then
TmpVal = SMALL
endif
- DumpBio(i,j,k) = TmpVal
+ DumpBio(i,j,k,1) = TmpVal
+ LimitCorr(i,j,k) = 2
! the positiveness is applied to
! the other components
@@ -170,6 +197,7 @@ subroutine wrt_chl_stat
MyCorr = bio%pquot(i,j,k,l)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,1)
MyCorr = MyCorr/LIM_THETA - bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
bio%phy(i,j,k,l,m) = max(0., MyCorr)
+ LimitCorr(i,j,k) = 3
endif
endif
@@ -182,6 +210,7 @@ subroutine wrt_chl_stat
MyCorr = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,2)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,2)
MyCorr = MyCorr*OPT_N_C - bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
bio%phy(i,j,k,l,m) = max(0., MyCorr)
+ LimitCorr(i,j,k) = 4
endif
endif
@@ -195,6 +224,7 @@ subroutine wrt_chl_stat
MyCorr = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,2)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,2)
MyCorr = MyCorr*OPT_P_C - bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
bio%phy(i,j,k,l,m) = max(0., MyCorr)
+ LimitCorr(i,j,k) = 5
endif
endif
@@ -208,16 +238,17 @@ subroutine wrt_chl_stat
MyCorr = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,2)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,2)
MyCorr = MyCorr*OPT_S_C - bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
bio%phy(i,j,k,l,m) = max(0., MyCorr)
+ LimitCorr(i,j,k) = 6
endif
endif
endif ! ApplyConditions
- DumpBio(i,j,k) = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,m)
+ DumpBio(i,j,k,1) = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,m)
endif
else
- DumpBio(i,j,k) = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
+ DumpBio(i,j,k,1) = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
endif
endif
@@ -225,7 +256,15 @@ subroutine wrt_chl_stat
enddo
enddo
- ierr = nf90mpi_put_var_all(ncid,idP,DumpBio,MyStart,MyCount)
+ do k=my_km+1,grd%km
+ do j=1,grd%jm
+ do i=1,grd%im
+ DumpBio(i,j,k,1) = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
+ enddo
+ enddo
+ enddo
+
+ ierr = nf90mpi_put_var_all(ncid,idP,DumpBio,MyStart_4d,MyCount_4d)
if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all '//DA_VarList(iVar), ierr)
ierr = nf90mpi_put_var_all(ncid,idTim,TimeArr,MyStartSingle,MyCountSingle)
@@ -243,6 +282,52 @@ subroutine wrt_chl_stat
enddo ! l
enddo ! m
+! File for post check DA
+! plus check variables
+ LimCorrfile = 'DA__FREQ_1/limcorr.'//ShortDate//'.nc'
+ ierr = nf90mpi_create(Var3DCommunicator, LimCorrfile, NF90_CLOBBER, MPI_INFO_NULL, ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('LimCorrfile ', ierr)
+
+ ierr = nf90mpi_def_dim(ncid,'x',global_im ,xid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim longitude ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'y' ,global_jm ,yid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim latitude ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'z' ,global_km, depid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim depth ', ierr)
+
+ LimVarName='lim_to_corr_FLAG'
+
+ ierr = nf90mpi_def_var(ncid, LimVarName, nf90_int, (/xid,yid,depid/), idL )
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', ierr)
+
+ ! ierr = nf90mpi_def_var(ncid, 'valtest', nf90_float, (/xid,yid/), iVar1 )
+ ! if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', ierr)
+
+ ! ierr = nf90mpi_def_var(ncid, 'initchl', nf90_float, (/xid,yid/), iVar2 )
+ ! if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', ierr)
+
+ ierr = nf90mpi_enddef(ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_enddef', ierr)
+
+ ierr = nf90mpi_put_var_all(ncid,idL,LimitCorr,MyStart,MyCount)
+ ! ierr = nf90mpi_put_var_all(ncid,idL,LimitCorr,MyStartLim,MyCountLim)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all ', ierr)
+
+ ! ierr =
+ ! nf90mpi_put_var_all(ncid,iVar1,ValuesToTest(:,:,1),MyStartLim,MyCountLim)
+ ! if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all ', ierr)
+
+ ! ierr =
+ ! nf90mpi_put_var_all(ncid,iVar2,bio%InitialChl(:,:,1),MyStartLim,MyCountLim)
+ ! if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all ', ierr)
+
+ ierr = nf90mpi_close(ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('LimCorrfile ', ierr)
+
+
+
+
DEALLOCATE(DumpBio, ValuesToTest, MyConditions)
-end subroutine wrt_chl_stat
\ No newline at end of file
+
+end subroutine wrt_chl_stat
diff --git a/wrt_dia.f90 b/wrt_dia.f90
index 4a98550..e33276b 100644
--- a/wrt_dia.f90
+++ b/wrt_dia.f90
@@ -47,7 +47,7 @@ subroutine wrt_dia
integer status
integer :: ncid,xid,yid,depid,idchl,idn3n,idn1p,ido2o
integer :: idvip,idmsk,eofid
- integer(kind=MPI_OFFSET_KIND) :: global_im, global_jm, global_km
+ integer(kind=MPI_OFFSET_KIND) :: global_im, global_jm, global_km, my_km
real(r4), allocatable, dimension(:,:,:) :: DumpMatrix
@@ -66,6 +66,10 @@ subroutine wrt_dia
global_im = GlobalRow
global_jm = GlobalCol
global_km = grd%km
+
+ my_km = grd%km
+ if(drv%multiv.eq.1) &
+ my_km = ros%kmchl
status = nf90mpi_def_dim(ncid,'depth' ,global_km, depid)
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim depth ', status)
@@ -74,24 +78,26 @@ subroutine wrt_dia
status = nf90mpi_def_dim(ncid,'longitude',global_im ,xid)
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim longitude ', status)
- if(drv%chl_assim .eq. 1) then
+ if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) then
status = nf90mpi_def_var(ncid,'chl', nf90_float, (/xid,yid,depid/), idchl )
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var chl', status)
status = nf90mpi_put_att(ncid,idchl , 'missing_value',1.e+20)
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_att', status)
endif
- if(drv%nut .eq. 1 .and. bio%n3n .eq. 1) then
+ if((drv%nut .eq. 1 .and. bio%n3n .eq. 1) .or. (drv%multiv .eq. 1)) then
status = nf90mpi_def_var(ncid,'n3n', nf90_float, (/xid,yid,depid/), idn3n )
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var n3n', status)
status = nf90mpi_put_att(ncid,idn3n , 'missing_value',1.e+20)
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_att', status)
endif
- if(drv%nut .eq. 1 .and. bio%n3n .eq. 1 .and. bio%updateN1p .eq. 1) then
+ if(bio%updateN1p .eq. 1) then
+ if((drv%nut .eq. 1 .and. bio%n3n .eq. 1) .or. (drv%multiv .eq. 1)) then
status = nf90mpi_def_var(ncid,'n1p', nf90_float, (/xid,yid,depid/), idn1p )
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var n1p', status)
status = nf90mpi_put_att(ncid,idn1p , 'missing_value',1.e+20)
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_att', status)
endif
+ endif
if(drv%nut .eq. 1 .and. bio%o2o .eq. 1) then
status = nf90mpi_def_var(ncid,'o2o', nf90_float, (/xid,yid,depid/), ido2o )
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var o2o', status)
@@ -102,8 +108,8 @@ subroutine wrt_dia
status = nf90mpi_enddef(ncid)
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', status)
- if(drv%chl_assim .eq. 1) then
- do k=1,grd%km
+ if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) then
+ do k=1,my_km
do j=1,grd%jm
do i=1,grd%im
if(grd%msk(i,j,k) .eq. 1) then
@@ -114,12 +120,23 @@ subroutine wrt_dia
enddo
enddo
enddo
+ do k=my_km+1,grd%km
+ do j=1,grd%jm
+ do i=1,grd%im
+ if(grd%msk(i,j,k) .eq. 1) then
+ DumpMatrix(i,j,k) = 0
+ else
+ DumpMatrix(i,j,k) = 1.e20
+ endif
+ enddo
+ enddo
+ enddo
status = nf90mpi_put_var_all(ncid,idchl,DumpMatrix,MyStart,MyCount)
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all chl', status)
endif
- if(drv%nut .eq. 1 .and. bio%n3n .eq. 1) then
+ if((drv%nut .eq. 1 .and. bio%n3n .eq. 1) .or. (drv%multiv .eq. 1)) then
do k=1,grd%km
do j=1,grd%jm
do i=1,grd%im
@@ -135,7 +152,8 @@ subroutine wrt_dia
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all n3n', status)
endif
- if(drv%nut .eq. 1 .and. bio%n3n .eq. 1 .and. bio%updateN1p .eq. 1) then
+ if (bio%updateN1p .eq. 1) then
+ if((drv%nut .eq. 1 .and. bio%n3n .eq. 1).or.(drv%multiv.eq.1)) then
do k=1,grd%km
do j=1,grd%jm
do i=1,grd%im
@@ -150,6 +168,7 @@ subroutine wrt_dia
status = nf90mpi_put_var_all(ncid,idn1p,DumpMatrix,MyStart,MyCount)
if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all n1p', status)
endif
+ endif
if(drv%nut .eq. 1 .and. bio%o2o .eq. 1) then
do k=1,grd%km
diff --git a/wrt_nut_stat.f90 b/wrt_nut_stat.f90
index ad0bf8d..5a78eec 100644
--- a/wrt_nut_stat.f90
+++ b/wrt_nut_stat.f90
@@ -17,31 +17,31 @@ subroutine wrt_nut_stat
INTEGER(kind=MPI_OFFSET_KIND) :: global_im, global_jm, global_km, MyTime
INTEGER(KIND=MPI_OFFSET_KIND) :: MyCountSingle(1), MyStartSingle(1)
- CHARACTER(LEN=37) :: BioRestart
- CHARACTER(LEN=39) :: BioRestartLong
+ CHARACTER(LEN=46) :: BioRestart
+ CHARACTER(LEN=47) :: BioRestartLong
CHARACTER(LEN=6) :: MyVarName
- LOGICAL, ALLOCATABLE :: MyConditions(:,:,:,:)
+ ! LOGICAL, ALLOCATABLE :: MyConditions(:,:,:,:)
- real(r8) :: TmpVal, MyCorr, MyRatio, SMALL
- real(r4), allocatable, dimension(:,:,:) :: DumpBio
+ real(r8) :: TmpVal, MyCorr, MyRatio!, SMALL
real(r4), allocatable, dimension(:,:,:,:) :: ValuesToTest
+
+ ! bug fix Intel 2018
+ real(r4), allocatable, dimension(:,:,:,:) :: DumpBio
+ integer(KIND=MPI_OFFSET_KIND) :: MyStart_4d(4), MyCount_4d(4)
+
real(r8) :: TimeArr(1)
- !real(r4) :: MAX_N_CHL, MAX_P_CHL, MAX_P_C, MAX_N_C
- !real(r4) :: OPT_N_C, OPT_P_C, OPT_S_C, LIM_THETA
-
- SMALL = 1.e-5
- ! MAX_N_CHL = 150. ! Derived from max chl:c=0.02 (BFMconsortium)
- ! MAX_P_CHL = 10.
- ! MAX_P_C = 7.86e-4*2 ! values from BFMconsortium parametrs document (P.Lazzari)
- ! OPT_P_C = 7.86e-4
- ! MAX_N_C = 1.26e-2*2 ! values from BFMconsortium parametrs document (P.Lazzari)
- ! OPT_N_C = 1.26e-2
- ! OPT_S_C = 0.01 ! values from BFMconsortium parametrs document (P.Lazzari)
- ! LIM_THETA = 0.01
+
+! SMALL = 1.e-5
+
+ MyStart_4d(1:3) = MyStart(:)
+ MyStart_4d(4) = 1
+ MyCount_4d(1:3) = MyCount(:)
+ MyCount_4d(4) = 1
+
- ALLOCATE(DumpBio(grd%im,grd%jm,grd%km)); DumpBio(:,:,:) = 1.e20
+ ALLOCATE(DumpBio(grd%im,grd%jm,grd%km,1)); DumpBio(:,:,:,:) = 1.e20
ALLOCATE(ValuesToTest(grd%im,grd%jm,grd%km,NNutVar)); ValuesToTest(:,:,:,:) = dble(0.)
- ALLOCATE(MyConditions(grd%im,grd%jm,grd%km,bio%nphy))
+ ! ALLOCATE(MyConditions(grd%im,grd%jm,grd%km,bio%nphy))
if(MyId .eq. 0) then
write(drv%dia,*) 'writing nut structure'
@@ -71,21 +71,8 @@ subroutine wrt_nut_stat
! if(bio%ApplyConditions) then
! !if(ValuesToTest(i,j,k) .gt. 10*bio%InitialChl(i,j,k)) then
- ! ! do m=1,bio%ncmp
- ! ! do l=1,bio%nphy
- ! ! bio%phy(i,j,k,1) = 9.*bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
- ! ! enddo
- ! ! enddo
-
! ! endif
- ! ! limitations in case of high nutrient contents
- ! do l=1,bio%nphy
- ! MyConditions(i,j,k,l) = bio%cquot(i,j,k,l,3) .gt. MAX_N_CHL
- ! MyConditions(i,j,k,l) = MyConditions(i,j,k,l) .or. (bio%cquot(i,j,k,l,4) .gt. MAX_P_CHL)
- ! MyConditions(i,j,k,l) = MyConditions(i,j,k,l) .or. (bio%cquot(i,j,k,l,3)/bio%cquot(i,j,k,l,2) .gt. (4*MAX_N_C))
- ! MyConditions(i,j,k,l) = MyConditions(i,j,k,l) .or. (bio%cquot(i,j,k,l,4)/bio%cquot(i,j,k,l,2) .gt. (4*MAX_P_C))
- ! enddo
! endif
endif
enddo
@@ -94,6 +81,7 @@ subroutine wrt_nut_stat
+
do l=1,NNutVar
iVar = NPhytoVar + l
@@ -102,8 +90,8 @@ subroutine wrt_nut_stat
write(*,*) "Warning: Reading a variable not in the DA_VarList!"
endif
- BioRestart = 'RESTARTS/RST.'//ShortDate//'.'//DA_VarList(iVar)//'.nc'
- BioRestartLong = 'RESTARTS/RST.'//DA_DATE//'.'//DA_VarList(iVar)//'.nc'
+ BioRestart = 'DA__FREQ_1/RST_after.'//ShortDate//'.'//DA_VarList(iVar)//'.nc'
+ BioRestartLong = 'DA__FREQ_1/RST_after.'//DA_DATE//'.'//DA_VarList(iVar)//'.nc'
if(drv%Verbose .eq. 1 .and. MyId .eq. 0) &
print*, "Writing Nut Restart ", BioRestart
@@ -146,76 +134,21 @@ subroutine wrt_nut_stat
! This correction must be the first
! condition applied (before apply corrections
! on the other components)
- TmpVal = 0.01*bio%InitialNut(i,j,k,l)
- if(TmpVal.gt.SMALL) then
- TmpVal = SMALL
- endif
- DumpBio(i,j,k) = TmpVal
+ TmpVal = 0.1*bio%InitialNut(i,j,k,l)
+ ! if(TmpVal.gt.SMALL) then
+ ! TmpVal = SMALL
+ ! endif
+ DumpBio(i,j,k,1) = TmpVal
else
- DumpBio(i,j,k) = ValuesToTest(i,j,k,l)
+ DumpBio(i,j,k,1) = ValuesToTest(i,j,k,l)
! if(bio%ApplyConditions) then
- ! if(bio%phy(i,j,k,l,m) .gt. 0 .and. MyConditions(i,j,k,l)) then
- ! bio%phy(i,j,k,l,m) = 0.
- ! endif
-
- ! ! limitation on Carbon corrections
- ! ! when chl/Carbon ratio is small
- ! if(m .eq. 2) then
- ! MyRatio = 1./bio%cquot(i,j,k,l,m)
- ! if(MyRatio .lt. LIM_THETA .and. bio%phy(i,j,k,l,m) .gt. 0) then
- ! MyCorr = bio%pquot(i,j,k,l)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,1)
- ! MyCorr = MyCorr/LIM_THETA - bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
- ! bio%phy(i,j,k,l,m) = max(0., MyCorr)
- ! endif
- ! endif
-
- ! ! limitation on Nitrogen corrections
- ! ! to the optimal N/C ratio
- ! if(m .eq. 3) then
- ! ! compute N/C fraction
- ! MyRatio = bio%cquot(i,j,k,l,m)/bio%cquot(i,j,k,l,2)
- ! if(MyRatio .gt. OPT_N_C .and. bio%phy(i,j,k,l,m) .gt. 0) then
- ! MyCorr = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,2)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,2)
- ! MyCorr = MyCorr*OPT_N_C - bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
- ! bio%phy(i,j,k,l,m) = max(0., MyCorr)
- ! endif
-
- ! endif
-
- ! ! limitation on Phosphorus corrections
- ! ! to the optimal P/C ratio
- ! if(m .eq. 4) then
- ! ! compute P/C fraction
- ! MyRatio = bio%cquot(i,j,k,l,m)/bio%cquot(i,j,k,l,2)
- ! if(MyRatio .gt. OPT_P_C .and. bio%phy(i,j,k,l,m) .gt. 0) then
- ! MyCorr = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,2)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,2)
- ! MyCorr = MyCorr*OPT_P_C - bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
- ! bio%phy(i,j,k,l,m) = max(0., MyCorr)
- ! endif
-
- ! endif
-
- ! ! limitation on Silicon corrections
- ! ! to the optimal Si/C ratio
- ! if(m .eq. 5) then
- ! ! compute Si/C fraction
- ! MyRatio = bio%cquot(i,j,k,l,m)/bio%cquot(i,j,k,l,2)
- ! if(MyRatio .gt. OPT_S_C .and. bio%phy(i,j,k,l,m) .gt. 0) then
- ! MyCorr = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,2)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,2)
- ! MyCorr = MyCorr*OPT_S_C - bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k)
- ! bio%phy(i,j,k,l,m) = max(0., MyCorr)
- ! endif
-
- ! endif
-
! endif ! ApplyConditions
- ! DumpBio(i,j,k) = bio%pquot(i,j,k,l)*bio%cquot(i,j,k,l,m)*bio%InitialChl(i,j,k) + bio%phy(i,j,k,l,m)
endif
else
- DumpBio(i,j,k) = bio%InitialNut(i,j,k,l)
+ DumpBio(i,j,k,1) = bio%InitialNut(i,j,k,l)
endif
endif
@@ -223,7 +156,7 @@ subroutine wrt_nut_stat
enddo
enddo
- ierr = nf90mpi_put_var_all(ncid,idP,DumpBio,MyStart,MyCount)
+ ierr = nf90mpi_put_var_all(ncid,idP,DumpBio,MyStart_4d,MyCount_4d)
if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all '//DA_VarList(iVar), ierr)
ierr = nf90mpi_put_var_all(ncid,idTim,TimeArr,MyStartSingle,MyCountSingle)
@@ -240,6 +173,7 @@ subroutine wrt_nut_stat
endif
enddo ! l
- DEALLOCATE(DumpBio, ValuesToTest, MyConditions)
+ DEALLOCATE(DumpBio, ValuesToTest)
+ ! DEALLOCATE(DumpBio, ValuesToTest, MyConditions)
-end subroutine wrt_nut_stat
\ No newline at end of file
+end subroutine wrt_nut_stat
diff --git a/wrt_upd_nut.f90 b/wrt_upd_nut.f90
new file mode 100644
index 0000000..42c7efb
--- /dev/null
+++ b/wrt_upd_nut.f90
@@ -0,0 +1,185 @@
+subroutine wrt_upd_nut
+
+ use set_knd
+ use grd_str
+ use drv_str
+ use mpi_str
+ use bio_str
+ use pnetcdf
+ use da_params
+
+ implicit none
+
+ INTEGER(i4) :: ncid, ierr, i, j, k, l, m
+ INTEGER(i4) :: idP, iVar
+ INTEGER(I4) :: xid,yid,depid,timeId, idTim
+ INTEGER :: system, SysErr
+
+ INTEGER(kind=MPI_OFFSET_KIND) :: global_im, global_jm, global_km, MyTime
+ INTEGER(KIND=MPI_OFFSET_KIND) :: MyCountSingle(1), MyStartSingle(1)
+ CHARACTER(LEN=45) :: BioRestart
+ CHARACTER(LEN=47) :: BioRestartLong
+ CHARACTER(LEN=6) :: MyVarName
+
+
+ real(r8) :: TmpVal, Myn3nUpdate, Myn1pUpdate, totchl, SMALL
+ real(r4), allocatable, dimension(:,:,:,:) :: ValuesToTest
+
+ ! bug fix Intel 2018
+ real(r4), allocatable, dimension(:,:,:,:) :: DumpBio
+ integer(KIND=MPI_OFFSET_KIND) :: MyStart_4d(4), MyCount_4d(4)
+
+ real(r8) :: TimeArr(1)
+
+ SMALL = 1.e-5
+
+ MyStart_4d(1:3) = MyStart(:)
+ MyStart_4d(4) = 1
+ MyCount_4d(1:3) = MyCount(:)
+ MyCount_4d(4) = 1
+
+
+ ALLOCATE(DumpBio(grd%im,grd%jm,grd%km,1)); DumpBio(:,:,:,:) = 1.e20
+ ALLOCATE(ValuesToTest(grd%im,grd%jm,grd%km,2)); ValuesToTest(:,:,:,:) = dble(0.)
+
+ if(MyId .eq. 0) then
+ write(drv%dia,*) 'writing nut update based on chl'
+ write(*,*) 'writing nut update based on chl'
+ endif
+
+ global_im = GlobalRow
+ global_jm = GlobalCol
+ global_km = grd%km
+ MyTime = 1
+
+ MyCountSingle(1) = 1
+ MyStartSingle(1) = 1
+ TimeArr(1) = DA_JulianDate
+
+
+ do k=1,grd%km
+ do j=1,grd%jm
+ do i=1,grd%im
+ totchl = 0.
+ do l=1,bio%nphy
+ totchl = totchl + bio%phy(i,j,k,l,1)
+ enddo
+ Myn3nUpdate = totchl*bio%covn3n_chl(i,j,k)
+ ValuesToTest(i,j,k,1) = bio%InitialNut(i,j,k,1) + Myn3nUpdate
+
+ Myn1pUpdate = totchl*bio%covn1p_chl(i,j,k)
+ ValuesToTest(i,j,k,2) = bio%InitialNut(i,j,k,2) + Myn1pUpdate
+
+ enddo
+ enddo
+ enddo
+
+ ! do k=1,grd%km
+ ! do j=1,grd%jm
+ ! do i=1,grd%im
+ ! totchl = 0.
+ ! do l=1,bio%nphy
+ ! totchl = totchl + bio%phy(i,j,k,l,1)
+ ! enddo
+ ! Myn3nUpdate = totchl*bio%covn3n_chl(i,j,k)
+ ! ValuesToTest(i,j,k,1) = bio%InitialNut(i,j,k,1) + Myn3nUpdate
+
+ ! ValuesToTest(i,j,k,2) = bio%InitialNut(i,j,k,2) + Myn3nUpdate*bio%covn3n_n1p(i,j,k)
+
+ ! enddo
+ ! enddo
+ ! enddo
+
+ do l=1,NNutVar
+ iVar = NPhytoVar + l
+
+ if(iVar .gt. NBioVar) then
+ if(MyId .eq. 0) &
+ write(*,*) "Warning: Reading a variable not in the DA_VarList!"
+ endif
+
+
+ BioRestart = 'DA__FREQ_1/RST_after.'//ShortDate//'.'//DA_VarList(iVar)//'.nc'
+ BioRestartLong = 'DA__FREQ_1/RST_after.'//DA_DATE//'.'//DA_VarList(iVar)//'.nc'
+
+ if(drv%Verbose .eq. 1 .and. MyId .eq. 0) &
+ print*, "Writing Nut Restart based on chl ", BioRestart
+
+ ierr = nf90mpi_create(Var3DCommunicator, BioRestart, NF90_CLOBBER, MPI_INFO_NULL, ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_create '//BioRestart, ierr)
+
+ ierr = nf90mpi_def_dim(ncid,'x',global_im ,xid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim longitude ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'y' ,global_jm ,yid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim latitude ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'z' ,global_km, depid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim depth ', ierr)
+ ierr = nf90mpi_def_dim(ncid,'time',MyTime ,timeId)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_dim time ', ierr)
+
+ MyVarName='TRN'//DA_VarList(iVar)
+
+ ierr = nf90mpi_def_var(ncid, MyVarName, nf90_float, (/xid,yid,depid,timeId/), idP )
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', ierr)
+
+ ierr = nf90mpi_def_var(ncid,'time' , nf90_double, (/timeid/) , idTim)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', ierr)
+ ierr = nf90mpi_put_att(ncid,idP , 'missing_value',1.e+20)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_att', ierr)
+
+ ierr = nf90mpi_enddef(ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_enddef'//DA_VarList(iVar), ierr)
+
+ do k=1,grd%km
+ do j=1,grd%jm
+ do i=1,grd%im
+
+ if(bio%InitialNut(i,j,k,l) .lt. 1.e20) then
+
+ if(grd%msk(i,j,k).eq.1) then
+
+ if(ValuesToTest(i,j,k,l) .lt. 0) then
+ ! Excluding negative concentrations
+ ! This correction must be the first
+ ! condition applied (before apply corrections
+ ! on the other components)
+ TmpVal = 0.01*bio%InitialNut(i,j,k,l)
+ if(TmpVal.gt.SMALL) then
+ TmpVal = SMALL
+ endif
+ DumpBio(i,j,k,1) = TmpVal
+
+ else
+ DumpBio(i,j,k,1) = ValuesToTest(i,j,k,l)
+ endif
+ else
+ DumpBio(i,j,k,1) = bio%InitialNut(i,j,k,l)
+ endif
+
+ endif
+ enddo
+ enddo
+ enddo
+
+ ierr = nf90mpi_put_var_all(ncid,idP,DumpBio,MyStart_4d,MyCount_4d)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all '//DA_VarList(iVar), ierr)
+
+ ierr = nf90mpi_put_var_all(ncid,idTim,TimeArr,MyStartSingle,MyCountSingle)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all '//DA_VarList(iVar), ierr)
+
+ ierr = nf90mpi_close(ncid)
+ if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_close '//BioRestart, ierr)
+
+ call MPI_Barrier(Var3DCommunicator, ierr)
+ ! only process 0 creates link to restart files
+ if(MyId .eq. 0) then
+ SysErr = system("ln -sf $PWD/"//BioRestart//" "//BioRestartLong)
+ if(SysErr /= 0) call MPI_Abort(MPI_COMM_WORLD, -1, SysErr)
+ endif
+ enddo ! l
+
+
+ DEALLOCATE(DumpBio, ValuesToTest)
+
+
+end subroutine wrt_upd_nut