From 97a3dc8148c9d067d8fd895f550b43d3ba357c77 Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Mon, 27 Oct 2025 16:34:25 +0100 Subject: [PATCH 01/13] Starting coupledDA and renaming oceanvar-->biovar --- CMakeLists.txt | 2 +- Makefile | 2 +- oceanvar.f90 => biovar.f90 | 4 ++-- def_nml_multi.f90 | 2 ++ main.f90 | 2 +- 5 files changed, 7 insertions(+), 5 deletions(-) rename oceanvar.f90 => biovar.f90 (99%) diff --git a/CMakeLists.txt b/CMakeLists.txt index c49e6d9..b2c37d7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,7 +26,7 @@ set(OBJS routines.f def_nml.f90 def_grd.f90 sav_itr.f90 ini_itr.f90 ini_bmd.f90 ini_nrm.f90 min_cfn.f90 costf.f90 cnv_ctv.f90 ver_hor.f90 rcfl_x.f90 rcfl_y.f90 veof.f90 obsop.f90 obs_sat.f90 resid.f90 res_inc.f90 rdrcorr.f90 mean_rdr.f90 obs_trd_ad.f90 obs_sat_ad.f90 veof_ad.f90 ver_hor_ad.f90 rcfl_x_ad.f90 rcfl_y_ad.f90 - cnv_ctv_ad.f90 cnv_inn.f90 wrt_dia.f90 filename_mod.f90 oceanvar.f90 main.f90) + cnv_ctv_ad.f90 cnv_inn.f90 wrt_dia.f90 filename_mod.f90 biovar.f90 main.f90) # SOURCES --END-- add_executable(var_3d ${KNDSTR} ${OBJSTR} ${OBJS}) diff --git a/Makefile b/Makefile index a9b6bdc..9ea3926 100644 --- a/Makefile +++ b/Makefile @@ -167,7 +167,7 @@ OBJS = \ readGrid.o\ def_cov.o\ tao_minimizer.o\ - oceanvar.o + biovar.o MAINEXE = main.o diff --git a/oceanvar.f90 b/biovar.f90 similarity index 99% rename from oceanvar.f90 rename to biovar.f90 index 2185806..1c40b7e 100644 --- a/oceanvar.f90 +++ b/biovar.f90 @@ -1,4 +1,4 @@ -subroutine oceanvar +subroutine biovar !--------------------------------------------------------------------------- ! ! @@ -146,4 +146,4 @@ subroutine oceanvar !----------------------------------------------------------------- if(MyId .eq. 0) close(drv%dia) -end subroutine oceanvar +end subroutine biovar diff --git a/def_nml_multi.f90 b/def_nml_multi.f90 index d5a0293..dfe4982 100644 --- a/def_nml_multi.f90 +++ b/def_nml_multi.f90 @@ -86,6 +86,7 @@ subroutine def_nml_multi 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,*) ' Copuled density/nutrient ass densnut = ', densnut 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 @@ -100,6 +101,7 @@ subroutine def_nml_multi drv%chl_upnut = chl_upnut drv%nut = nut drv%multiv = multiv + drv%densnut = densnut bio%nphy = nphyto sat%dep = chl_dep bio%ncmp = ncmp diff --git a/main.f90 b/main.f90 index af4dded..7d2fa98 100644 --- a/main.f90 +++ b/main.f90 @@ -12,7 +12,7 @@ program ocean_var ! in case of standalone usage call var3d_mpi_init -call oceanvar +call biovar ! finalizing the MPI environment call clean_da_params From 143cabcaf9c9774f9691434d2a4cc3a28edffeff Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Fri, 31 Oct 2025 12:12:17 +0100 Subject: [PATCH 02/13] Some PHY parts removed and adding densnut parameter --- CMakeLists.txt | 7 ------- Makefile | 37 ------------------------------------- def_nml_multi.f90 | 8 ++++---- drv_str.f90 | 1 + obs_str.f90 | 2 +- 5 files changed, 6 insertions(+), 49 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b2c37d7..a024865 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,13 +13,6 @@ set(KNDSTR set_knd.f90) set(OBJSTR drv_str.f90 cns_str.f90 obs_str.f90 grd_str.f90 bmd_str.f90 eof_str.f90 ctl_str.f90) -set(PHYSOBS bar_mod.f90 bar_mod_ad.f90 div_dmp.f90 div_dmp_ad.f90 get_byg.f90 get_byg_ad.f90 - get_obs_arg.f90 get_obs_gld.f90 get_obs_xbt.f90 get_obs_sla.f90 - get_obs_vdr.f90 get_obs_gvl.f90 get_obs_tra.f90 get_obs_trd.f90 - obs_arg.f90 obs_xbt.f90 obs_gld.f90 obs_vdr.f90 obs_gvl.f90 obs_tra.f90 obs_trd.f90 obs_sla.f90 - obsop_ad.f90 obs_sla_ad.f90 obs_arg_ad.f90 obs_xbt_ad.f90 obs_gld_ad.f90 obs_vdr_ad.f90 obs_gvl_ad.f90 obs_tra_ad.f90 - mod_trj_tl.f90 mod_trj_ad.f90 get_vel_ad.f90 get_vel.f90 - invrt.f90 invrt_ad.f90 ) set(OBJS routines.f def_nml.f90 def_grd.f90 sav_itr.f90 ini_itr.f90 ini_bmd.f90 rdgrds.f90 rdeofs.f90 netcdf_err.f90 get_obs.f90get_obs_sat.f90 int_par.f90 obs_vec.f90 def_cov.f90 ini_cfn.f90 diff --git a/Makefile b/Makefile index 9ea3926..1c7b792 100644 --- a/Makefile +++ b/Makefile @@ -62,43 +62,6 @@ OBJSTR = \ bio_str.o\ mpi_str.o -PHYSOBS = \ - get_obs_sla.o\ - get_obs_arg.o\ - get_obs_xbt.o\ - get_obs_gld.o\ - get_obs_vdr.o\ - get_obs_gvl.o\ - get_obs_tra.o\ - get_obs_trd.o\ - obs_sla.o\ - obs_arg.o\ - obs_xbt.o\ - obs_gld.o\ - obs_vdr.o\ - obs_gvl.o\ - obs_tra.o\ - obs_trd.o\ - obs_sla_ad.o\ - obs_arg_ad.o\ - obs_xbt_ad.o\ - obs_gld_ad.o\ - obs_vdr_ad.o\ - obs_gvl_ad.o\ - obs_tra_ad.o\ - obs_trd_ad.o\ - bar_mod.o\ - get_vel.o\ - div_dmp.o\ - get_byg.o\ - bar_mod_ad.o\ - get_vel_ad.o\ - div_dmp_ad.o\ - get_byg_ad.o\ - mod_trj_ad.o\ - mod_trj_tl.o\ - invrt.o\ - invrt_ad.o OBJS = \ def_nml.o\ diff --git a/def_nml_multi.f90 b/def_nml_multi.f90 index dfe4982..265a6cd 100644 --- a/def_nml_multi.f90 +++ b/def_nml_multi.f90 @@ -47,14 +47,14 @@ subroutine def_nml_multi LOGICAL :: ApplyConditions INTEGER(i4) :: chl_assim, chl_upnut, nut, multiv, N3n, O2o, updateN1p - INTEGER(i4) :: nphyto, uniformL, anisL + INTEGER(i4) :: nphyto, uniformL, anisL, densnut 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 + NAMELIST /params/ sat_obs, argo, densnut, uniformL, anisL ! ------------------------------------------------------------------- @@ -86,7 +86,6 @@ subroutine def_nml_multi 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,*) ' Copuled density/nutrient ass densnut = ', densnut 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 @@ -101,7 +100,6 @@ subroutine def_nml_multi drv%chl_upnut = chl_upnut drv%nut = nut drv%multiv = multiv - drv%densnut = densnut bio%nphy = nphyto sat%dep = chl_dep bio%ncmp = ncmp @@ -121,6 +119,7 @@ subroutine def_nml_multi 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,*) ' Copuled density/nutrient ass densnut = ', densnut write(drv%dia,*) ' Set uniform correlation radius uniformL = ', uniformL write(drv%dia,*) ' Set anisotropy on corr radius anisL = ', anisL write(drv%dia,*) '------------------------------------------------------------' @@ -134,6 +133,7 @@ subroutine def_nml_multi drv%sat_obs = sat_obs drv%argo_obs = argo + drv%densnut = densnut drv%uniformL = uniformL drv%anisL = anisL diff --git a/drv_str.f90 b/drv_str.f90 index 3e73d84..40826f5 100644 --- a/drv_str.f90 +++ b/drv_str.f90 @@ -41,6 +41,7 @@ MODULE drv_str INTEGER(i4) :: MyCounter ! Number of iteration done by Tao solver 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) :: densnut ! Flag for the update of nutrients based on density increments 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) diff --git a/obs_str.f90 b/obs_str.f90 index 9d5eee7..c522393 100644 --- a/obs_str.f90 +++ b/obs_str.f90 @@ -92,7 +92,7 @@ MODULE obs_str TYPE (arg_t) :: arg ! --- - ! Observational vector for Chlorophyll + ! Observational vector for SAT chl TYPE chl_t INTEGER(i8) :: no ! Number of all observations From ccfe9f284061260b74c7702b02184df833dbc0df Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Sat, 1 Nov 2025 23:37:04 +0100 Subject: [PATCH 03/13] Preparing new files and strategy for DA Not working!!! --- cns_str.f90 | 4 - cnv_inn.f90 | 16 +- filename_mod.f90 | 4 +- get_phyincr_arg.f90 | 410 ++++++++++++++++++++++++++++++++++++++++++++ ver_hor_densnut.f90 | 251 +++++++++++++++++++++++++++ 5 files changed, 675 insertions(+), 10 deletions(-) create mode 100644 get_phyincr_arg.f90 create mode 100644 ver_hor_densnut.f90 diff --git a/cns_str.f90 b/cns_str.f90 index 2d87c46..3e215e3 100644 --- a/cns_str.f90 +++ b/cns_str.f90 @@ -39,13 +39,11 @@ MODULE cns_str INTEGER(i4) :: ntr ! No. of iterations (half of) REAL(r8) :: dx ! Grid resolution (m) REAL(r8) :: L ! Correlation radius -!laura REAL(r8),POINTER :: Lxyz(:,:,:)!Correlation radius from file in km REAL(r8),POINTER :: L_x(:,:,:)!Correlation radius from file in km REAL(r8),POINTER :: L_y(:,:,:)!Correlation radius from file in km REAL(r8),POINTER :: rtx(:,:)!Correlation radius from file in km REAL(r8),POINTER :: rty(:,:)!Correlation radius from file in km -!laura REAL(r8) :: E ! Norm REAL(r8) :: alp ! Filter weight INTEGER(i4) :: ntb ! Number of points in the table @@ -53,9 +51,7 @@ MODULE cns_str REAL(r8) :: dsmx ! Maximum distance REAL(r8) :: dsl ! Table increment REAL(r8), POINTER :: al(:) ! Filter weights in the table -!laura REAL(r8), POINTER :: sc(:,:) ! Filter scaling factors in the table -!laura REAL(r8) :: scl ! Scaling factor REAL(r8) :: efc ! Scaling factor for extended points INTEGER(i4) :: kstp ! Step for extended points diff --git a/cnv_inn.f90 b/cnv_inn.f90 index ca686f0..fc91d59 100644 --- a/cnv_inn.f90 +++ b/cnv_inn.f90 @@ -49,17 +49,25 @@ subroutine cnv_inn 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%densnut .eq. 1) then + call ver_hor_densnut(grd%n3n,'N') + else + call ver_hor_nut(grd%n3n, grd%n3n_ad,'N') + endif endif if(bio%O2o .eq. 1) then - call ver_hor_nut(grd%o2o, grd%o2o_ad,'O') + 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') + if (drv%densnut .eq. 1) then + call ver_hor_densnut(grd%n3n,'N') + else + call ver_hor_nut(grd%n3n, grd%n3n_ad,'N') + endif endif ! --- diff --git a/filename_mod.f90 b/filename_mod.f90 index dcf6e14..a91a9ce 100644 --- a/filename_mod.f90 +++ b/filename_mod.f90 @@ -14,10 +14,10 @@ MODULE FILENAMES character (LEN=1024) :: EIV_FILE != 'eiv.nc' character (LEN=1024) :: OBS_FILE != 'obs_1.dat' character (LEN=1024) :: GRID_FILE != 'grid1.nc' -!laura character (LEN=1024) :: RCORR_FILE != 'chl_rad_corr.nc' character (LEN=1024) :: ARGO_FILE != 'argo_mis.dat' character (LEN=1024) :: ANIS_FILE != 'gradsal.nc' +character (LEN=1024) :: INCR_FILE != 'density_increments.nc' CONTAINS @@ -38,10 +38,10 @@ SUBROUTINE SETFILENAMES EIV_FILE = 'eiv.nc' OBS_FILE = 'obs_1.dat' ! 'obs_'//fgrd//'.dat' GRID_FILE = 'grid1.nc'! 'grid'//cgrd//'.nc' -!laura RCORR_FILE = 'chl_rad_corr.nc' ARGO_FILE = 'arg_mis.dat' ANIS_FILE = 'gradsal.nc' + INCR_FILE = 'density_increments.dat' END SUBROUTINE SETFILENAMES diff --git a/get_phyincr_arg.f90 b/get_phyincr_arg.f90 new file mode 100644 index 0000000..7a45d6d --- /dev/null +++ b/get_phyincr_arg.f90 @@ -0,0 +1,410 @@ +subroutine get_obs_arg + + !--------------------------------------------------------------------------- + ! ! + ! 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 . ! + ! ! + !--------------------------------------------------------------------------- + + !----------------------------------------------------------------------- + ! ! + ! Load ARGO observations ! + ! ! + ! Version 1: S.Dobricic 2006 ! + !----------------------------------------------------------------------- + + use set_knd + use drv_str + use grd_str + use obs_str + use mpi_str + use filenames + use bio_str + + implicit none + + 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 + + arg%no = 0 + arg%nc = 0 + + + ! --- + ! Allocate memory for observations + if(MyId .eq. 0) then + open(511,file=trim(ARGO_FILE)) + read(511,'(I4)') GlobalArgNum + write(drv%dia,*)'Number of ARGO observations: ', GlobalArgNum + endif + + call MPI_Bcast(GlobalArgNum, 1, MPI_INT, 0, Var3DCommunicator, ierr) + + if(GlobalArgNum .eq. 0)then + if(MyId .eq. 0) & + close(511) + return + endif + + ALLOCATE( TmpFlc(GlobalArgNum), TmpPar(GlobalArgNum)) + 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 + ! process 0 reads all the argo observations + do k=1,GlobalArgNum + read (511,*) & + TmpFlc(k), TmpPar(k), & + TmpLon(k), TmpLat(k), & + TmpDpt(k), TmpTim(k), & + TmpRes(k), TmpErr(k), & + ! TmpStd(k), TmpIno(k) + TmpIno(k) + end do + close (511) + endif + + call MPI_Bcast(TmpFlc, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpPar, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpLon, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpLat, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpDpt, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpTim, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpRes, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpErr, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpIno, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + + ! Counting the number of observations that falls in the domain + Counter = 0 + 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.(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 + enddo + + if(drv%Verbose .eq. 1) & + print*, "MyId", MyId, "has",Counter,"ARGO observations" + + arg%no = Counter + + ALLOCATE ( arg%ino(arg%no), arg%flg(arg%no), arg%flc(arg%no), arg%par(arg%no)) + ALLOCATE ( arg%lon(arg%no), arg%lat(arg%no), arg%dpt(arg%no), arg%tim(arg%no)) + 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)) + ALLOCATE ( arg%pq5(arg%no), arg%pq6(arg%no), arg%pq7(arg%no), arg%pq8(arg%no)) + + Counter = 0 + 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. 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) + arg%par(Counter) = TmpPar(k) + arg%lon(Counter) = TmpLon(k) + arg%lat(Counter) = TmpLat(k) + 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 + + enddo + + ! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST + ! arg%res(:) = 1 + ! arg%err(:) = 1.d-2 + + + ! --- + ! Initialise quality flag + arg%flg(:) = 1 + arg%rb(:) = 0 + + ! --- +! Vertical interpolation parameters + do k = 1,arg%no + if(arg%flg(k).eq.1)then + arg%kb(k) = 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 + arg%kb(k) = kk + arg%rb(k) = (arg%dpt(k) - grd%dep(kk)) / (grd%dep(kk+1) - grd%dep(kk)) + else if ( arg%dpt(k).ge.0 .and. arg%dpt(k).lt.grd%dep(1)) then + arg%kb(k) = 1 + endif + enddo + endif + enddo + + + ! --- + ! Count good observations + arg%nc = 0 + do k=1,arg%no + if(arg%flg(k).eq.1)then + arg%nc = arg%nc + 1 + else + arg%res(k) = 0. + arg%inc(k) = 0. + arg%pq1(k) = 0. + arg%pq2(k) = 0. + arg%pq3(k) = 0. + arg%pq4(k) = 0. + arg%pq5(k) = 0. + arg%pq6(k) = 0. + arg%pq7(k) = 0. + arg%pq8(k) = 0. + endif + enddo + arg%flc(:) = arg%flg(:) + + DEALLOCATE( TmpFlc, TmpPar) + DEALLOCATE( TmpLon, TmpLat) + DEALLOCATE( TmpDpt, TmpTim) + DEALLOCATE( TmpRes, TmpErr) +! DEALLOCATE( TMPStd, TmpIno) + DEALLOCATE( TmpIno) + +end subroutine get_obs_arg + + + +subroutine int_par_arg + + !----------------------------------------------------------------------- + ! ! + ! Get interpolation parameters for a grid ! + ! ! + ! Version 1: S.Dobricic 2006 ! + !----------------------------------------------------------------------- + + 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, kind, kk + INTEGER(i4) :: i1, j1, k1, idep + REAL(r8) :: p1, q1, r1 + REAL(r8) :: msk4, div_x, div_y + LOGICAL :: ins + + ins(i,i1) = i.ge.1 .and. i.le.i1 + + if(arg%no.gt.0) then + + arg%flc(:) = arg%flg(:) + + ! --- + ! Horizontal interpolation parameters + do k = 1,arg%no + do j=1,grd%jm-1 + do i=1,grd%im-1 + if( grd%lat(i,j).le.arg%lat(k) .and. grd%lat(i,j+1).gt.arg%lat(k) .and. & + grd%lon(i,j).le.arg%lon(k) .and. grd%lon(i+1,j).gt.arg%lon(k) ) then + j1 = j + i1 = i + q1 = j1 + (arg%lat(k) - grd%lat(i,j)) / (grd%lat(i,j+1) - grd%lat(i,j)) + p1 = i1 + (arg%lon(k) - grd%lon(i,j)) / (grd%lon(i+1,j) - grd%lon(i,j)) + else if( i .eq. grd%im-1 .and. grd%lat(i,j).le.arg%lat(k) .and. grd%lat(i,j+1).gt.arg%lat(k) .and. & + grd%lon(grd%im,j).le.arg%lon(k) .and. grd%NextLongitude.gt.arg%lon(k) ) then + j1 = j + i1 = grd%im + q1 = j1 + (arg%lat(k) - grd%lat(i,j)) / (grd%lat(i,j+1) - grd%lat(i,j)) + p1 = i1 + (arg%lon(k) - grd%lon(grd%im,j)) / (grd%NextLongitude - grd%lon(grd%im,j)) + endif + enddo + enddo + + ! q1 = (arg%lat(k) - grd%lat(1,1)) / grd%dlt + 1.0 + ! j1 = int(q1) + ! p1 = (arg%lon(k) - grd%lon(1,1)) / grd%dln + 1.0 + ! i1 = int(p1) + if(ins(j1,grd%jm) .and. ins(i1,grd%im)) then + arg%ib(k) = i1 + arg%jb(k) = j1 + arg%pb(k) = (p1-i1) + arg%qb(k) = (q1-j1) + else + arg%flc(k) = 0 + endif + enddo + + ! --- + ! Undefine masked for multigrid + do k = 1,arg%no + if(arg%flc(k).eq.1)then + i1 = arg%ib(k) + j1 = arg%jb(k) + idep = arg%kb(k)+1 + msk4 = grd%global_msk(GlobalRowOffset+i1,j1,idep) + grd%global_msk(GlobalRowOffset+i1+1,j1,idep) + & + grd%global_msk(GlobalRowOffset+i1,j1+1,idep) + grd%global_msk(GlobalRowOffset+i1+1,j1+1,idep) + if(msk4.lt.1.) arg%flc(k) = 0 + endif + enddo + + ! --- + ! Horizontal interpolation parameters for each masked grid + do k = 1,arg%no + 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) + j1=arg%jb(k) + q1=arg%qb(k) + + + k1=arg%kb(k) + div_y = (1.-q1) * max(grd%global_msk(GlobalRowOffset+i1,j1 ,k1),grd%global_msk(GlobalRowOffset+i1+1,j1 ,k1)) & + + q1 * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) + div_x = (1.-p1) * grd%global_msk(GlobalRowOffset+i1 ,j1,k1) + p1 * grd%global_msk(GlobalRowOffset+i1+1,j1,k1) + arg%pq1(k) = grd%global_msk(GlobalRowOffset+i1,j1,k1) & + * max(grd%global_msk(GlobalRowOffset+i1,j1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1,k1)) & + * (1.-p1) * (1.-q1) & + /( div_x * div_y + 1.e-16 ) + arg%pq2(k) = grd%global_msk(GlobalRowOffset+i1+1,j1,k1) & + * max(grd%global_msk(GlobalRowOffset+i1,j1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1,k1)) & + * p1 * (1.-q1) & + /( div_x * div_y + 1.e-16 ) + div_x = (1.-p1) * grd%global_msk(GlobalRowOffset+i1 ,j1+1,k1) + p1 * grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) + arg%pq3(k) = grd%global_msk(GlobalRowOffset+i1,j1+1,k1) & + * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) & + * (1.-p1) * q1 & + /( div_x * div_y + 1.e-16 ) + arg%pq4(k) = grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) & + * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) & + * p1 * q1 & + /( div_x * div_y + 1.e-16 ) + + k1=arg%kb(k) + 1 + div_y = (1.-q1) * max(grd%global_msk(GlobalRowOffset+i1,j1 ,k1),grd%global_msk(GlobalRowOffset+i1+1,j1 ,k1)) & + + q1 * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) + div_x = (1.-p1) * grd%global_msk(GlobalRowOffset+i1 ,j1,k1) + p1 * grd%global_msk(GlobalRowOffset+i1+1,j1,k1) + arg%pq5(k) = grd%global_msk(GlobalRowOffset+i1,j1,k1) & + * max(grd%global_msk(GlobalRowOffset+i1,j1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1,k1)) & + * (1.-p1) * (1.-q1) & + /( div_x * div_y + 1.e-16 ) + arg%pq6(k) = grd%global_msk(GlobalRowOffset+i1+1,j1,k1) & + * max(grd%global_msk(GlobalRowOffset+i1,j1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1,k1)) & + * p1 * (1.-q1) & + /( div_x * div_y + 1.e-16 ) + div_x = (1.-p1) * grd%global_msk(GlobalRowOffset+i1 ,j1+1,k1) + p1 * grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) + arg%pq7(k) = grd%global_msk(GlobalRowOffset+i1,j1+1,k1) & + * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) & + * (1.-p1) * q1 & + /( div_x * div_y + 1.e-16 ) + arg%pq8(k) = grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) & + * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) & + * p1 * q1 & + /( div_x * div_y + 1.e-16 ) + + r1=arg%rb(k) + arg%pq1(k) = (1.-r1) * arg%pq1(k) + arg%pq2(k) = (1.-r1) * arg%pq2(k) + arg%pq3(k) = (1.-r1) * arg%pq3(k) + arg%pq4(k) = (1.-r1) * arg%pq4(k) + arg%pq5(k) = r1 * arg%pq5(k) + arg%pq6(k) = r1 * arg%pq6(k) + arg%pq7(k) = r1 * arg%pq7(k) + arg%pq8(k) = r1 * arg%pq8(k) + + if(arg%pq1(k) .lt. 1.E-16) arg%pq1(k) = dble(0) + if(arg%pq2(k) .lt. 1.E-16) arg%pq2(k) = dble(0) + if(arg%pq3(k) .lt. 1.E-16) arg%pq3(k) = dble(0) + if(arg%pq4(k) .lt. 1.E-16) arg%pq4(k) = dble(0) + if(arg%pq5(k) .lt. 1.E-16) arg%pq5(k) = dble(0) + if(arg%pq6(k) .lt. 1.E-16) arg%pq6(k) = dble(0) + if(arg%pq7(k) .lt. 1.E-16) arg%pq7(k) = dble(0) + if(arg%pq8(k) .lt. 1.E-16) arg%pq8(k) = dble(0) + + + 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 + arg%nc = 0 + do k=1,arg%no + if(arg%flc(k).eq.1)then + arg%nc = arg%nc + 1 + endif + enddo + + endif + + arg%nc_global = 0 + call MPI_Allreduce(arg%nc, arg%nc_global, 1, MPI_INT, MPI_SUM, Var3DCommunicator, ierr) + + if(MyId .eq. 0) then + write(drv%dia,*)'Real number of ARGO observations: ',arg%nc_global + print*,'Good argo observations: ',arg%nc_global + end if + + DEALLOCATE ( arg%ino, arg%flg) + DEALLOCATE ( arg%lon, arg%lat, arg%dpt, arg%tim) + DEALLOCATE ( arg%pb, arg%qb, arg%rb) + +end subroutine int_par_arg diff --git a/ver_hor_densnut.f90 b/ver_hor_densnut.f90 new file mode 100644 index 0000000..be99ae2 --- /dev/null +++ b/ver_hor_densnut.f90 @@ -0,0 +1,251 @@ +subroutine ver_hor_nut(NutArray, NutArrayAd, Var) + + !--------------------------------------------------------------------------- + ! ! + ! 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 . ! + ! ! + !--------------------------------------------------------------------------- + + !----------------------------------------------------------------------- + ! ! + ! Apply horizontal filter ! + ! ! + ! Version 1: S.Dobricic 2006 ! + ! Version 2: S.Dobricic 2007 ! + ! Symmetric calculation in presence of coastal boundaries ! + ! eta_ad, tem_ad, and sal_ad are here temporary arrays ! + ! Version 3: A. Teruzzi 2013 ! + ! Attenuation of correction near the cost where d<200m ! + !----------------------------------------------------------------------- + + + use set_knd + use grd_str + use eof_str + use cns_str + use drv_str + use obs_str + use mpi_str + + implicit none + + INTEGER(i4) :: i,j,k, ione, l + INTEGER :: jp, SurfaceIndex, TmpOffset, LinearIndex + INTEGER(i4) :: iProc, ierr + type(DoubleGrid), allocatable, dimension(:,:,:) :: SendBuf3D + type(DoubleGrid), allocatable, dimension(:) :: RecBuf1D(:) + REAL(r8), allocatable, dimension(:,:,:) :: DefBufChl, DefBufChlAd + REAL(r8) :: NutArray(grd%im,grd%jm,grd%km), NutArrayAd(grd%im,grd%jm,grd%km) + CHARACTER :: Var + + ione = 1 + + ! --- + ! Vertical EOFs + call veof_nut(NutArray, Var) + !return + ! goto 103 !No Vh + + ! --- + ! Load temporary arrays + do k=1,grd%km + NutArrayAd(:,:,k) = NutArray(:,:,k) + enddo + + !********** APPLY RECURSIVE FILTERS ********** ! + ! --- + ! Transpose calculation in the presense of coastal boundaries + + ! --- + ! y direction + ! --- + ! Scale by the scaling factor + do k=1,grd%km + NutArrayAd(:,:,k) = NutArrayAd(:,:,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, NutArrayAd, grd%jnx, grd%jmx) + + ! --- + ! 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)) + + do k=1,grd%km + do j=1,grd%jm + do i=1,grd%im + SendBuf3D(k,i,j)%chl = NutArray(i,j,k) + end do + end do + end do + do k=1,grd%km + do j=1,grd%jm + do i=1,grd%im + SendBuf3D(k,i,j)%chl_ad = NutArrayAd(i,j,k) + end do + end do + end do + + call MPI_Alltoallv(SendBuf3D, SendCountX3D, SendDisplX3D, MyPair, & + RecBuf1D, RecCountX3D, RecDisplX3D, MyPair, Var3DCommunicator, ierr) + + SurfaceIndex = localCol*grd%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 + DefBufChl(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl + end do + end do + + end do + 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 + DefBufChlAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl_ad + end do + end do + end do + end do + + ! --- + ! Scale by the scaling factor + do k=1,grd%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) + + else + ! --- + ! Scale by the scaling factor + do k=1,grd%km + NutArrayAd(:,:,k) = NutArrayAd(:,:,k) * grd%scx(:,:,k) + enddo + + call rcfl_x_ad( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, NutArrayAd, grd%inx, grd%imx) + + end if + + + + ! --- + ! 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) + + do k=1,grd%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)) + + do k=1,grd%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 j=1,localCol + do i=1,GlobalRow + SendBuf3D(k,j,i)%chl_ad = DefBufChlAd(i,j,k) + end do + end do + end do + + call MPI_Alltoallv(SendBuf3D, RecCountX3D, RecDisplX3D, MyPair, & + RecBuf1D, SendCountX3D, SendDisplX3D, MyPair, Var3DCommunicator, ierr) + + SurfaceIndex = grd%im*grd%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 + NutArray(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl + end do + end do + end do + 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 + NutArrayAd(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl_ad + end do + end do + end do + end do + + DEALLOCATE(SendBuf3D, RecBuf1D, DefBufChl, DefBufChlAd) + + else ! NumProcI .eq. 1 + + call rcfl_x( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, NutArray, grd%inx, grd%imx) + + do k=1,grd%km + NutArray(:,:,k) = NutArray(:,:,k) * grd%scx(:,:,k) + enddo + + end if + + ! --- + ! y direction + ! Apply recursive filter in y direction + call rcfl_y( localRow, GlobalCol, grd%km, grd%jmax, grd%aey, grd%bey, NutArray, grd%jnx, grd%jmx) + + ! --- + ! Scale by the scaling factor + do k=1,grd%km + NutArray(:,:,k) = NutArray(:,:,k) * grd%scy(:,:,k) + enddo + + ! --- + ! Average + do k=1,grd%km + NutArray(:,:,k) = (NutArray(:,:,k) + NutArrayAd(:,:,k) ) * 0.5 + enddo + + ! --- + ! Scale for boundaries + do k=1,grd%km + NutArray(:,:,k) = NutArray(:,:,k) * grd%msk(:,:,k) + enddo + + ! 103 continue + +end subroutine ver_hor_nut From fc4b5dc65fd2ae0c9e020ea7b40d573c0736d3de Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Wed, 5 Nov 2025 09:46:42 +0100 Subject: [PATCH 04/13] First version of assimilation of density increment profiles and some DEALLOCATE added --- Makefile | 6 +- biovar.f90 | 4 + clean_mem.f90 | 13 + cnv_inn.f90 | 20 +- def_cov.f90 | 6 +- dnc_str.f90 | 81 +++++ drv_str.f90 | 1 + get_phyincr_arg.f90 => get_densincr_arg.f90 | 361 +++++++++----------- get_obs.f90 | 4 + int_par.f90 | 8 +- rdeofs_multi.f90 | 2 +- readGrid.f90 | 3 + sav_itr.f90 | 14 + v_densnut.f90 | 91 +++++ ver_hor_densnut.f90 | 15 +- wrt_chl_stat.f90 | 2 +- wrt_dia.f90 | 62 +++- wrt_nut_stat.f90 | 4 +- 18 files changed, 468 insertions(+), 229 deletions(-) create mode 100644 dnc_str.f90 rename get_phyincr_arg.f90 => get_densincr_arg.f90 (50%) create mode 100644 v_densnut.f90 diff --git a/Makefile b/Makefile index 1c7b792..adaf875 100644 --- a/Makefile +++ b/Makefile @@ -60,7 +60,8 @@ OBJSTR = \ ctl_str.o\ rcfl_mod.o\ bio_str.o\ - mpi_str.o + mpi_str.o\ + dnc_str.o\ OBJS = \ @@ -78,16 +79,19 @@ OBJS = \ get_obs.o\ get_obs_arg.o\ get_obs_sat.o\ + get_densincr_arg.o\ int_par.o\ obs_vec.o\ ini_cfn.o\ cnv_ctv.o\ ver_hor_chl.o\ ver_hor_nut.o\ + ver_hor_densnut.o\ rcfl_x.o\ rcfl_y.o\ veof_chl.o\ veof_nut.o\ + v_densnut.o\ obsop.o\ obs_arg.o\ resid.o\ diff --git a/biovar.f90 b/biovar.f90 index 1c40b7e..6cf781e 100644 --- a/biovar.f90 +++ b/biovar.f90 @@ -137,6 +137,10 @@ subroutine biovar endif endif + if (drv%dnc .eq. 1) then + call wrt_nut_stat + endif + call sav_itr if(MyId .eq. 0) write(drv%dia,*) 'out of sav_itr ' diff --git a/clean_mem.f90 b/clean_mem.f90 index a05255b..355735d 100644 --- a/clean_mem.f90 +++ b/clean_mem.f90 @@ -39,6 +39,7 @@ subroutine clean_mem use cns_str use rcfl use mpi_str + use dnc_str implicit none @@ -68,6 +69,18 @@ subroutine clean_mem DEALLOCATE ( sat%dzr) endif + ! density increments + if(drv%dnc .eq. 1) then + DEALLOCATE ( dnc%flc) + DEALLOCATE ( dnc%inc) + DEALLOCATE ( dnc%corr) + DEALLOCATE ( dnc%err) + DEALLOCATE ( dnc%std) + DEALLOCATE ( dnc%ib, dnc%jb, dnc%kb) + DEALLOCATE ( dnc%pq1, dnc%pq2, dnc%pq3, dnc%pq4) + DEALLOCATE ( dnc%pq5, dnc%pq6, dnc%pq7, dnc%pq8) + endif + ! Constants structure DEALLOCATE ( rcf%al) DEALLOCATE ( rcf%sc) diff --git a/cnv_inn.f90 b/cnv_inn.f90 index fc91d59..83dcffc 100644 --- a/cnv_inn.f90 +++ b/cnv_inn.f90 @@ -49,27 +49,23 @@ subroutine cnv_inn endif if(drv%nut .eq. 1) then if(bio%N3n .eq. 1) then - if (drv%densnut .eq. 1) then - call ver_hor_densnut(grd%n3n,'N') - else - call ver_hor_nut(grd%n3n, grd%n3n_ad,'N') - endif + 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') + call ver_hor_nut(grd%o2o, grd%o2o_ad,'O') endif endif endif - + if (drv%multiv .eq. 1) then call ver_hor_chl - if (drv%densnut .eq. 1) then - call ver_hor_densnut(grd%n3n,'N') - else - call ver_hor_nut(grd%n3n, grd%n3n_ad,'N') - endif + call ver_hor_nut(grd%n3n, grd%n3n_ad,'N') endif + if (drv%densnut .eq. 1) then + grd%n3n(:,:,:) = 0.0 + call ver_hor_densnut(grd%n3n) + endif ! --- ! Apply biological repartition of the chlorophyll if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) & diff --git a/def_cov.f90 b/def_cov.f90 index bf81477..0fd2d81 100644 --- a/def_cov.f90 +++ b/def_cov.f90 @@ -545,6 +545,10 @@ subroutine def_cov call readNutStat if(bio%updateN1p.eq.1) & call readNutCov - endif + endif + if(drv%dnc .eq. 1) then + call readNutStat + call readNutCov + endif end subroutine def_cov diff --git a/dnc_str.f90 b/dnc_str.f90 new file mode 100644 index 0000000..da66cd6 --- /dev/null +++ b/dnc_str.f90 @@ -0,0 +1,81 @@ +MODULE dnc_str + + !--------------------------------------------------------------------------- + ! ! + ! 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 . ! + ! ! + !--------------------------------------------------------------------------- + + !----------------------------------------------------------------------- + ! ! + ! Density increment vectors ! + ! ! + ! Version 1: S.Dobricic 2006 ! + !----------------------------------------------------------------------- + + use set_knd + + implicit none + + public + + ! --- + ! Density increment vector for ARGO floats + TYPE dnc_t + + INTEGER(i8) :: no ! Number of all observations + INTEGER(i8) :: nc ! Number of good observations + INTEGER(i8) :: k ! Density increment index + ! REAL(r8) :: dep ! Minimum depth for observations + ! INTEGER(i8) :: kdp ! Model level corresponding to dep + ! INTEGER(i8), POINTER :: ino(:) ! Float number + ! INTEGER(i8), POINTER :: par(:) ! Parameter flag (0-chl, 1-N3n, 2-O2o) + INTEGER(i8), POINTER :: flg(:) ! Quality flag + INTEGER(i8), POINTER :: flc(:) ! Temporary flag for multigrid + REAL(r8), POINTER :: lon(:) ! Longitute + REAL(r8), POINTER :: lat(:) ! Latitude + REAL(r8), POINTER :: dpt(:) ! Depth + ! REAL(r8), POINTER :: tim(:) ! Time + REAL(r8), POINTER :: inc(:) ! Increments + REAL(r8), POINTER :: corr(:) ! Correlations + REAL(r8), POINTER :: err(:) ! Nitrate std (error) + REAL(r8), POINTER :: std(:) ! Density std + ! 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 + INTEGER(i8), POINTER :: jb(:) ! j index of the nearest south point + REAL(r8) , POINTER :: qb(:) ! distance from the nearest south point + INTEGER(i8), POINTER :: kb(:) ! k index of the nearest point below + REAL(r8) , POINTER :: rb(:) ! distance from the nearest point below + REAL(r8) , POINTER :: pq1(:) ! Interpolation parameter for masked grids + REAL(r8) , POINTER :: pq2(:) ! Interpolation parameter for masked grids + REAL(r8) , POINTER :: pq3(:) ! Interpolation parameter for masked grids + REAL(r8) , POINTER :: pq4(:) ! Interpolation parameter for masked grids + REAL(r8) , POINTER :: pq5(:) ! Interpolation parameter for masked grids + REAL(r8) , POINTER :: pq6(:) ! Interpolation parameter for masked grids + REAL(r8) , POINTER :: pq7(:) ! Interpolation parameter for masked grids + REAL(r8) , POINTER :: pq8(:) ! Interpolation parameter for masked grids + + INTEGER(i4) :: nc_global ! Number of global good observations + + END TYPE dnc_t + + TYPE (dnc_t) :: dnc + + +END MODULE dnc_str diff --git a/drv_str.f90 b/drv_str.f90 index 40826f5..46d5492 100644 --- a/drv_str.f90 +++ b/drv_str.f90 @@ -49,6 +49,7 @@ MODULE drv_str INTEGER(i4) :: Verbose ! Flag for printing verbose output INTEGER(i4) :: nut ! Flag for nutrient assimilation INTEGER(i4) :: multiv ! Flag for multivariate assimilation + INTEGER(i4) :: dnc ! Flag for update of nitrate based on density increments END TYPE drv_t diff --git a/get_phyincr_arg.f90 b/get_densincr_arg.f90 similarity index 50% rename from get_phyincr_arg.f90 rename to get_densincr_arg.f90 index 7a45d6d..0ea307c 100644 --- a/get_phyincr_arg.f90 +++ b/get_densincr_arg.f90 @@ -1,4 +1,4 @@ -subroutine get_obs_arg +subroutine get_densincr_arg !--------------------------------------------------------------------------- ! ! @@ -23,7 +23,7 @@ subroutine get_obs_arg !----------------------------------------------------------------------- ! ! - ! Load ARGO observations ! + ! Load Density increments ! ! ! ! Version 1: S.Dobricic 2006 ! !----------------------------------------------------------------------- @@ -31,144 +31,137 @@ subroutine get_obs_arg use set_knd use drv_str use grd_str - use obs_str + use dnc_str +! use obs_str use mpi_str use filenames - use bio_str +! use bio_str implicit none 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 + REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpFlc, TmpLon, TmpLat + REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpDpt, TmpErr, TmpStd + REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpInc, TmpCorr + INTEGER(i4) :: GlobalDncNum, Counter, ierr character(len=1024) :: filename - arg%no = 0 - arg%nc = 0 + dnc%no = 0 + dnc%nc = 0 ! --- ! Allocate memory for observations if(MyId .eq. 0) then - open(511,file=trim(ARGO_FILE)) - read(511,'(I4)') GlobalArgNum - write(drv%dia,*)'Number of ARGO observations: ', GlobalArgNum + open(511,file=trim(INCR_FILE)) + read(511,'(I4)') GlobalDncNum + write(drv%dia,*)'Number of density increments: ', GlobalDncNum endif - call MPI_Bcast(GlobalArgNum, 1, MPI_INT, 0, Var3DCommunicator, ierr) + call MPI_Bcast(GlobalDncNum, 1, MPI_INT, 0, Var3DCommunicator, ierr) - if(GlobalArgNum .eq. 0)then + if(GlobalDncNum .eq. 0)then if(MyId .eq. 0) & close(511) return endif - ALLOCATE( TmpFlc(GlobalArgNum), TmpPar(GlobalArgNum)) - ALLOCATE( TmpLon(GlobalArgNum), TmpLat(GlobalArgNum)) - ALLOCATE( TmpDpt(GlobalArgNum), TmpTim(GlobalArgNum)) - ALLOCATE( TmpRes(GlobalArgNum), TmpErr(GlobalArgNum)) -! ALLOCATE( TmpStd(GlobalArgNum), TmpIno(GlobalArgNum)) - ALLOCATE( TmpIno(GlobalArgNum)) + ALLOCATE( TmpFlc(GlobalDncNum)) !, TmpPar(GlobalDncNum)) + ALLOCATE( TmpLon(GlobalDncNum), TmpLat(GlobalDncNum)) + ALLOCATE( TmpDpt(GlobalDncNum))!, TmpTim(GlobalDncNum)) + ALLOCATE( TmpInc(GlobalDncNum), TmpCorr(GlobalDncNum)) + ALLOCATE( TmpErr(GlobalDncNum), TmpStd(GlobalDncNum)) if(MyId .eq. 0) then - ! process 0 reads all the argo observations - do k=1,GlobalArgNum + ! process 0 reads all the density increments + do k=1,GlobalDncNum read (511,*) & - TmpFlc(k), TmpPar(k), & + TmpFlc(k), & !TmpPar(k), & TmpLon(k), TmpLat(k), & - TmpDpt(k), TmpTim(k), & - TmpRes(k), TmpErr(k), & - ! TmpStd(k), TmpIno(k) - TmpIno(k) + TmpDpt(k), &!TmpTim(k), & + TmpInc(k), TmpCorr(k), & + TmpErr(k), TmpStd(k) end do close (511) endif - call MPI_Bcast(TmpFlc, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpPar, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpLon, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpLat, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpDpt, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpTim, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpRes, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpErr, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpIno, GlobalArgNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpFlc, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) +! call MPI_Bcast(TmpPar, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpLon, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpLat, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpDpt, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) +! call MPI_Bcast(TmpTim, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpInc, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpCorr, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpErr, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + call MPI_Bcast(TmpStd, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) ! Counting the number of observations that falls in the domain Counter = 0 - do k=1,GlobalArgNum + do k=1,GlobalDncNum 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.(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 + if(drv%dnc.eq.1) then Counter = Counter + 1 endif endif enddo if(drv%Verbose .eq. 1) & - print*, "MyId", MyId, "has",Counter,"ARGO observations" + print*, "MyId", MyId, "has",Counter,"Density increments" - arg%no = Counter + dnc%no = Counter - ALLOCATE ( arg%ino(arg%no), arg%flg(arg%no), arg%flc(arg%no), arg%par(arg%no)) - ALLOCATE ( arg%lon(arg%no), arg%lat(arg%no), arg%dpt(arg%no), arg%tim(arg%no)) - 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)) - ALLOCATE ( arg%pq5(arg%no), arg%pq6(arg%no), arg%pq7(arg%no), arg%pq8(arg%no)) + ALLOCATE ( dnc%flg(dnc%no), dnc%flc(dnc%no)) !, dnc%par(dnc%no)) + ALLOCATE ( dnc%lon(dnc%no), dnc%lat(dnc%no), dnc%dpt(dnc%no)) !, dnc%tim(dnc%no)) + ALLOCATE ( dnc%inc(dnc%no)) + ALLOCATE ( dnc%corr(dnc%no)) + ALLOCATE ( dnc%err(dnc%no)) + ALLOCATE ( dnc%std(dnc%no)) + ALLOCATE ( dnc%ib(dnc%no), dnc%jb(dnc%no), dnc%kb(dnc%no)) + ALLOCATE ( dnc%pb(dnc%no), dnc%qb(dnc%no), dnc%rb(dnc%no)) + ALLOCATE ( dnc%pq1(dnc%no), dnc%pq2(dnc%no), dnc%pq3(dnc%no), dnc%pq4(dnc%no)) + ALLOCATE ( dnc%pq5(dnc%no), dnc%pq6(dnc%no), dnc%pq7(dnc%no), dnc%pq8(dnc%no)) Counter = 0 - do k=1,GlobalArgNum + do k=1,GlobalDncNum 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. 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 + if(drv%dnc.eq.1) then Counter = Counter + 1 - arg%flc(Counter) = TmpFlc(k) - arg%par(Counter) = TmpPar(k) - arg%lon(Counter) = TmpLon(k) - arg%lat(Counter) = TmpLat(k) - 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) + dnc%flc(Counter) = TmpFlc(k) + ! dnc%par(Counter) = TmpPar(k) + dnc%lon(Counter) = TmpLon(k) + dnc%lat(Counter) = TmpLat(k) + dnc%dpt(Counter) = TmpDpt(k) + dnc%inc(Counter) = TmpInc(k) + dnc%corr(Counter) = TmpCorr(k) + dnc%err(Counter) = TmpErr(k) + dnc%std(Counter) = TmpStd(k) + ! dnc%ino(Counter) = TmpIno(k) endif endif enddo - ! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST - ! arg%res(:) = 1 - ! arg%err(:) = 1.d-2 - ! --- ! Initialise quality flag - arg%flg(:) = 1 - arg%rb(:) = 0 + dnc%flg(:) = 1 + dnc%rb(:) = 0 ! --- ! Vertical interpolation parameters - do k = 1,arg%no - if(arg%flg(k).eq.1)then - arg%kb(k) = grd%km-1 + do k = 1,dnc%no + if(dnc%flg(k).eq.1)then + dnc%kb(k) = 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 - arg%kb(k) = kk - arg%rb(k) = (arg%dpt(k) - grd%dep(kk)) / (grd%dep(kk+1) - grd%dep(kk)) - else if ( arg%dpt(k).ge.0 .and. arg%dpt(k).lt.grd%dep(1)) then - arg%kb(k) = 1 + if( dnc%dpt(k).ge.grd%dep(kk) .and. dnc%dpt(k).lt.grd%dep(kk+1) ) then + dnc%kb(k) = kk + dnc%rb(k) = (dnc%dpt(k) - grd%dep(kk)) / (grd%dep(kk+1) - grd%dep(kk)) + else if ( dnc%dpt(k).ge.0 .and. dnc%dpt(k).lt.grd%dep(1)) then + dnc%kb(k) = 1 endif enddo endif @@ -177,37 +170,38 @@ subroutine get_obs_arg ! --- ! Count good observations - arg%nc = 0 - do k=1,arg%no - if(arg%flg(k).eq.1)then - arg%nc = arg%nc + 1 + dnc%nc = 0 + do k=1,dnc%no + if(dnc%flg(k).eq.1)then + dnc%nc = dnc%nc + 1 else - arg%res(k) = 0. - arg%inc(k) = 0. - arg%pq1(k) = 0. - arg%pq2(k) = 0. - arg%pq3(k) = 0. - arg%pq4(k) = 0. - arg%pq5(k) = 0. - arg%pq6(k) = 0. - arg%pq7(k) = 0. - arg%pq8(k) = 0. + dnc%inc(k) = 0. + dnc%pq1(k) = 0. + dnc%pq2(k) = 0. + dnc%pq3(k) = 0. + dnc%pq4(k) = 0. + dnc%pq5(k) = 0. + dnc%pq6(k) = 0. + dnc%pq7(k) = 0. + dnc%pq8(k) = 0. endif enddo - arg%flc(:) = arg%flg(:) + dnc%flc(:) = dnc%flg(:) - DEALLOCATE( TmpFlc, TmpPar) + DEALLOCATE( TmpFlc)!, TmpPar) DEALLOCATE( TmpLon, TmpLat) - DEALLOCATE( TmpDpt, TmpTim) - DEALLOCATE( TmpRes, TmpErr) -! DEALLOCATE( TMPStd, TmpIno) - DEALLOCATE( TmpIno) + DEALLOCATE( TmpDpt)!, TmpTim) + DEALLOCATE( TmpInc)!, TmpTim) + DEALLOCATE( TmpCorr)!, TmpTim) + DEALLOCATE( TmpErr) + DEALLOCATE( TMPStd) +! DEALLOCATE( TmpIno) -end subroutine get_obs_arg +end subroutine get_densincr_arg -subroutine int_par_arg +subroutine int_par_dnc !----------------------------------------------------------------------- ! ! @@ -219,8 +213,9 @@ subroutine int_par_arg use set_knd use drv_str use grd_str - use eof_str - use obs_str +! use eof_str +! use obs_str + use dnc_str use mpi_str implicit none @@ -233,178 +228,160 @@ subroutine int_par_arg ins(i,i1) = i.ge.1 .and. i.le.i1 - if(arg%no.gt.0) then + if(dnc%no.gt.0) then - arg%flc(:) = arg%flg(:) + dnc%flc(:) = dnc%flg(:) ! --- ! Horizontal interpolation parameters - do k = 1,arg%no + do k = 1,dnc%no do j=1,grd%jm-1 do i=1,grd%im-1 - if( grd%lat(i,j).le.arg%lat(k) .and. grd%lat(i,j+1).gt.arg%lat(k) .and. & - grd%lon(i,j).le.arg%lon(k) .and. grd%lon(i+1,j).gt.arg%lon(k) ) then + if( grd%lat(i,j).le.dnc%lat(k) .and. grd%lat(i,j+1).gt.dnc%lat(k) .and. & + grd%lon(i,j).le.dnc%lon(k) .and. grd%lon(i+1,j).gt.dnc%lon(k) ) then j1 = j i1 = i - q1 = j1 + (arg%lat(k) - grd%lat(i,j)) / (grd%lat(i,j+1) - grd%lat(i,j)) - p1 = i1 + (arg%lon(k) - grd%lon(i,j)) / (grd%lon(i+1,j) - grd%lon(i,j)) - else if( i .eq. grd%im-1 .and. grd%lat(i,j).le.arg%lat(k) .and. grd%lat(i,j+1).gt.arg%lat(k) .and. & - grd%lon(grd%im,j).le.arg%lon(k) .and. grd%NextLongitude.gt.arg%lon(k) ) then + q1 = j1 + (dnc%lat(k) - grd%lat(i,j)) / (grd%lat(i,j+1) - grd%lat(i,j)) + p1 = i1 + (dnc%lon(k) - grd%lon(i,j)) / (grd%lon(i+1,j) - grd%lon(i,j)) + else if( i .eq. grd%im-1 .and. grd%lat(i,j).le.dnc%lat(k) .and. grd%lat(i,j+1).gt.dnc%lat(k) .and. & + grd%lon(grd%im,j).le.dnc%lon(k) .and. grd%NextLongitude.gt.dnc%lon(k) ) then j1 = j i1 = grd%im - q1 = j1 + (arg%lat(k) - grd%lat(i,j)) / (grd%lat(i,j+1) - grd%lat(i,j)) - p1 = i1 + (arg%lon(k) - grd%lon(grd%im,j)) / (grd%NextLongitude - grd%lon(grd%im,j)) + q1 = j1 + (dnc%lat(k) - grd%lat(i,j)) / (grd%lat(i,j+1) - grd%lat(i,j)) + p1 = i1 + (dnc%lon(k) - grd%lon(grd%im,j)) / (grd%NextLongitude - grd%lon(grd%im,j)) endif enddo enddo - ! q1 = (arg%lat(k) - grd%lat(1,1)) / grd%dlt + 1.0 + ! q1 = (dnc%lat(k) - grd%lat(1,1)) / grd%dlt + 1.0 ! j1 = int(q1) - ! p1 = (arg%lon(k) - grd%lon(1,1)) / grd%dln + 1.0 + ! p1 = (dnc%lon(k) - grd%lon(1,1)) / grd%dln + 1.0 ! i1 = int(p1) if(ins(j1,grd%jm) .and. ins(i1,grd%im)) then - arg%ib(k) = i1 - arg%jb(k) = j1 - arg%pb(k) = (p1-i1) - arg%qb(k) = (q1-j1) + dnc%ib(k) = i1 + dnc%jb(k) = j1 + dnc%pb(k) = (p1-i1) + dnc%qb(k) = (q1-j1) else - arg%flc(k) = 0 + dnc%flc(k) = 0 endif enddo ! --- ! Undefine masked for multigrid - do k = 1,arg%no - if(arg%flc(k).eq.1)then - i1 = arg%ib(k) - j1 = arg%jb(k) - idep = arg%kb(k)+1 + do k = 1,dnc%no + if(dnc%flc(k).eq.1)then + i1 = dnc%ib(k) + j1 = dnc%jb(k) + idep = dnc%kb(k)+1 msk4 = grd%global_msk(GlobalRowOffset+i1,j1,idep) + grd%global_msk(GlobalRowOffset+i1+1,j1,idep) + & grd%global_msk(GlobalRowOffset+i1,j1+1,idep) + grd%global_msk(GlobalRowOffset+i1+1,j1+1,idep) - if(msk4.lt.1.) arg%flc(k) = 0 + if(msk4.lt.1.) dnc%flc(k) = 0 endif enddo ! --- ! Horizontal interpolation parameters for each masked grid - do k = 1,arg%no - if(arg%flg(k) .eq. 1) then - ! if(arg%flg(k) .eq. 1) then ! to verify that it works also in this case + do k = 1,dnc%no + if(dnc%flg(k) .eq. 1) then + ! if(dnc%flg(k) .eq. 1) then ! to verify that it works also in this case - i1=arg%ib(k) - p1=arg%pb(k) - j1=arg%jb(k) - q1=arg%qb(k) + i1=dnc%ib(k) + p1=dnc%pb(k) + j1=dnc%jb(k) + q1=dnc%qb(k) - k1=arg%kb(k) + k1=dnc%kb(k) div_y = (1.-q1) * max(grd%global_msk(GlobalRowOffset+i1,j1 ,k1),grd%global_msk(GlobalRowOffset+i1+1,j1 ,k1)) & + q1 * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) div_x = (1.-p1) * grd%global_msk(GlobalRowOffset+i1 ,j1,k1) + p1 * grd%global_msk(GlobalRowOffset+i1+1,j1,k1) - arg%pq1(k) = grd%global_msk(GlobalRowOffset+i1,j1,k1) & + dnc%pq1(k) = grd%global_msk(GlobalRowOffset+i1,j1,k1) & * max(grd%global_msk(GlobalRowOffset+i1,j1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1,k1)) & * (1.-p1) * (1.-q1) & /( div_x * div_y + 1.e-16 ) - arg%pq2(k) = grd%global_msk(GlobalRowOffset+i1+1,j1,k1) & + dnc%pq2(k) = grd%global_msk(GlobalRowOffset+i1+1,j1,k1) & * max(grd%global_msk(GlobalRowOffset+i1,j1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1,k1)) & * p1 * (1.-q1) & /( div_x * div_y + 1.e-16 ) div_x = (1.-p1) * grd%global_msk(GlobalRowOffset+i1 ,j1+1,k1) + p1 * grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) - arg%pq3(k) = grd%global_msk(GlobalRowOffset+i1,j1+1,k1) & + dnc%pq3(k) = grd%global_msk(GlobalRowOffset+i1,j1+1,k1) & * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) & * (1.-p1) * q1 & /( div_x * div_y + 1.e-16 ) - arg%pq4(k) = grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) & + dnc%pq4(k) = grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) & * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) & * p1 * q1 & /( div_x * div_y + 1.e-16 ) - k1=arg%kb(k) + 1 + k1=dnc%kb(k) + 1 div_y = (1.-q1) * max(grd%global_msk(GlobalRowOffset+i1,j1 ,k1),grd%global_msk(GlobalRowOffset+i1+1,j1 ,k1)) & + q1 * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) div_x = (1.-p1) * grd%global_msk(GlobalRowOffset+i1 ,j1,k1) + p1 * grd%global_msk(GlobalRowOffset+i1+1,j1,k1) - arg%pq5(k) = grd%global_msk(GlobalRowOffset+i1,j1,k1) & + dnc%pq5(k) = grd%global_msk(GlobalRowOffset+i1,j1,k1) & * max(grd%global_msk(GlobalRowOffset+i1,j1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1,k1)) & * (1.-p1) * (1.-q1) & /( div_x * div_y + 1.e-16 ) - arg%pq6(k) = grd%global_msk(GlobalRowOffset+i1+1,j1,k1) & + dnc%pq6(k) = grd%global_msk(GlobalRowOffset+i1+1,j1,k1) & * max(grd%global_msk(GlobalRowOffset+i1,j1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1,k1)) & * p1 * (1.-q1) & /( div_x * div_y + 1.e-16 ) div_x = (1.-p1) * grd%global_msk(GlobalRowOffset+i1 ,j1+1,k1) + p1 * grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) - arg%pq7(k) = grd%global_msk(GlobalRowOffset+i1,j1+1,k1) & + dnc%pq7(k) = grd%global_msk(GlobalRowOffset+i1,j1+1,k1) & * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) & * (1.-p1) * q1 & /( div_x * div_y + 1.e-16 ) - arg%pq8(k) = grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) & + dnc%pq8(k) = grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1) & * max(grd%global_msk(GlobalRowOffset+i1,j1+1,k1),grd%global_msk(GlobalRowOffset+i1+1,j1+1,k1)) & * p1 * q1 & /( div_x * div_y + 1.e-16 ) - r1=arg%rb(k) - arg%pq1(k) = (1.-r1) * arg%pq1(k) - arg%pq2(k) = (1.-r1) * arg%pq2(k) - arg%pq3(k) = (1.-r1) * arg%pq3(k) - arg%pq4(k) = (1.-r1) * arg%pq4(k) - arg%pq5(k) = r1 * arg%pq5(k) - arg%pq6(k) = r1 * arg%pq6(k) - arg%pq7(k) = r1 * arg%pq7(k) - arg%pq8(k) = r1 * arg%pq8(k) + r1=dnc%rb(k) + dnc%pq1(k) = (1.-r1) * dnc%pq1(k) + dnc%pq2(k) = (1.-r1) * dnc%pq2(k) + dnc%pq3(k) = (1.-r1) * dnc%pq3(k) + dnc%pq4(k) = (1.-r1) * dnc%pq4(k) + dnc%pq5(k) = r1 * dnc%pq5(k) + dnc%pq6(k) = r1 * dnc%pq6(k) + dnc%pq7(k) = r1 * dnc%pq7(k) + dnc%pq8(k) = r1 * dnc%pq8(k) - if(arg%pq1(k) .lt. 1.E-16) arg%pq1(k) = dble(0) - if(arg%pq2(k) .lt. 1.E-16) arg%pq2(k) = dble(0) - if(arg%pq3(k) .lt. 1.E-16) arg%pq3(k) = dble(0) - if(arg%pq4(k) .lt. 1.E-16) arg%pq4(k) = dble(0) - if(arg%pq5(k) .lt. 1.E-16) arg%pq5(k) = dble(0) - if(arg%pq6(k) .lt. 1.E-16) arg%pq6(k) = dble(0) - if(arg%pq7(k) .lt. 1.E-16) arg%pq7(k) = dble(0) - if(arg%pq8(k) .lt. 1.E-16) arg%pq8(k) = dble(0) + if(dnc%pq1(k) .lt. 1.E-16) dnc%pq1(k) = dble(0) + if(dnc%pq2(k) .lt. 1.E-16) dnc%pq2(k) = dble(0) + if(dnc%pq3(k) .lt. 1.E-16) dnc%pq3(k) = dble(0) + if(dnc%pq4(k) .lt. 1.E-16) dnc%pq4(k) = dble(0) + if(dnc%pq5(k) .lt. 1.E-16) dnc%pq5(k) = dble(0) + if(dnc%pq6(k) .lt. 1.E-16) dnc%pq6(k) = dble(0) + if(dnc%pq7(k) .lt. 1.E-16) dnc%pq7(k) = dble(0) + if(dnc%pq8(k) .lt. 1.E-16) dnc%pq8(k) = dble(0) 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 - arg%nc = 0 - do k=1,arg%no - if(arg%flc(k).eq.1)then - arg%nc = arg%nc + 1 + dnc%nc = 0 + do k=1,dnc%no + if(dnc%flc(k).eq.1)then + dnc%nc = dnc%nc + 1 endif enddo endif - arg%nc_global = 0 - call MPI_Allreduce(arg%nc, arg%nc_global, 1, MPI_INT, MPI_SUM, Var3DCommunicator, ierr) + dnc%nc_global = 0 + call MPI_Allreduce(dnc%nc, dnc%nc_global, 1, MPI_INT, MPI_SUM, Var3DCommunicator, ierr) if(MyId .eq. 0) then - write(drv%dia,*)'Real number of ARGO observations: ',arg%nc_global - print*,'Good argo observations: ',arg%nc_global + write(drv%dia,*)'Real number of denisty increments: ',dnc%nc_global + print*,'Good denisty increments: ',dnc%nc_global end if - DEALLOCATE ( arg%ino, arg%flg) - DEALLOCATE ( arg%lon, arg%lat, arg%dpt, arg%tim) - DEALLOCATE ( arg%pb, arg%qb, arg%rb) + DEALLOCATE ( dnc%flg) + DEALLOCATE ( dnc%lon, dnc%lat, dnc%dpt) !, dnc%tim) + DEALLOCATE ( dnc%pb, dnc%qb, dnc%rb) -end subroutine int_par_arg +end subroutine int_par_dnc diff --git a/get_obs.f90 b/get_obs.f90 index ed087a3..b506641 100644 --- a/get_obs.f90 +++ b/get_obs.f90 @@ -51,4 +51,8 @@ subroutine get_obs if(drv%sat_obs .eq. 1) & call get_obs_sat + ! Load density increments for Physical DA + if(drv%densnut .eq. 1) & + call get_densincr_arg + end subroutine get_obs diff --git a/int_par.f90 b/int_par.f90 index c6a533c..685d215 100644 --- a/int_par.f90 +++ b/int_par.f90 @@ -44,5 +44,9 @@ subroutine int_par ! Load observations of chlorophyll if(drv%sat_obs .eq. 1) & call int_par_chl - -end subroutine int_par + + ! Load density increments for Physical DA + if(drv%densnut .eq. 1) & + call int_par_dnc + + end subroutine int_par diff --git a/rdeofs_multi.f90 b/rdeofs_multi.f90 index 2dfcbc7..11bb373 100644 --- a/rdeofs_multi.f90 +++ b/rdeofs_multi.f90 @@ -205,7 +205,7 @@ subroutine rdeofs_multi stat = nf90mpi_close(ncid) if (stat /= nf90_noerr) call handle_err("nf90mpi_close", stat) - DEALLOCATE(x3, x2, std_chl, std_n3n) + DEALLOCATE(x3, x2, x1, std_chl, std_n3n) end subroutine rdeofs_multi diff --git a/readGrid.f90 b/readGrid.f90 index 28ef216..25ff1f3 100644 --- a/readGrid.f90 +++ b/readGrid.f90 @@ -162,6 +162,9 @@ subroutine readGrid 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 + if (drv%dnc .eq. 1) then + ALLOCATE (grd%n3n(grd%im,grd%jm,grd%km) ) ; grd%n3n = huge(grd%n3n(1,1,1)) + endif 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)) diff --git a/sav_itr.f90 b/sav_itr.f90 index 6df2aff..f5c8e56 100644 --- a/sav_itr.f90 +++ b/sav_itr.f90 @@ -41,6 +41,7 @@ subroutine sav_itr use mpi_str use bio_str use da_params + use dnc_str implicit none @@ -91,6 +92,10 @@ subroutine sav_itr DEALLOCATE( grd%n3n) DEALLOCATE( grd%n3n_ad) endif + + if(drv%dnc .eq. 1) then + DEALLOCATE( grd%n3n) + endif ! Observational vector DEALLOCATE( obs%inc, obs%amo, obs%res) @@ -111,6 +116,10 @@ subroutine sav_itr endif endif + if(drv%multiv.eq.1) then + DEALLOCATE( ros%evc_multi, ros%eva_multi) + endif + ! Control structure DEALLOCATE( ctl%x_c, ctl%g_c) @@ -147,6 +156,11 @@ subroutine sav_itr if(bio%updateN1p.eq.1) DEALLOCATE( bio%covn3n_n1p) endif + if(drv%dnc .eq. 1) then + DEALLOCATE( bio%InitialNut) + DEALLOCATE( bio%covn3n_n1p) + endif + DEALLOCATE(SurfaceWaterPoints) DEALLOCATE( a_rcx) diff --git a/v_densnut.f90 b/v_densnut.f90 new file mode 100644 index 0000000..e33a4c5 --- /dev/null +++ b/v_densnut.f90 @@ -0,0 +1,91 @@ +subroutine v_densnut(NutArray) + + + !--------------------------------------------------------------------------- + ! ! + ! 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 . ! + ! ! + !--------------------------------------------------------------------------- + + !----------------------------------------------------------------------- + ! ! + ! Apply density nutrient increments and unterpolare aon vertical grid ! + ! ! + ! Version 1: A. Teruzzi 2025 ! + !----------------------------------------------------------------------- + + + use set_knd + use grd_str + ! use eof_str + ! use obs_str + use mpi_str + ! use filenames + use drv_str + use bio_str + use dnc_str + + implicit none + + INTEGER(i4) :: i, j, k, kk, my_km + REAL(r8), DIMENSION(grd%jm,grd%km) :: slicevar + REAL(r8), DIMENSION(grd%im+1,grd%jm,grd%km) :: N3nExt + REAL(r8) :: NutArray(grd%im,grd%jm,grd%km) + REAL(8) :: nit_dinc + + my_km = grd%km + + ! call EXTEND_2D( grd%n3n_ad, grd%km, N3nExtended_3d ) + + dnc%k = 0 + N3nExt(:,:,:) = 0.0 + + do kk = 1,dnc%no + + i=dnc%ib(kk) + j=dnc%jb(kk) + k=dnc%kb(kk) + + if(dnc%flc(kk).eq.1) then + + dnc%k = dnc%k + 1 + nit_dinc = dnc%inc(dnc%k) * dnc%corr(dnc%k) * dnc%err(dnc%k) / dnc%std(dnc%k) + + N3nExt(i ,j ,k ) = N3nExt(i ,j ,k ) + dnc%pq1(kk) * nit_dinc + N3nExt(i+1,j ,k ) = N3nExt(i+1,j ,k ) + dnc%pq2(kk) * nit_dinc + N3nExt(i ,j+1,k ) = N3nExt(i ,j+1,k ) + dnc%pq3(kk) * nit_dinc + N3nExt(i+1,j+1,k ) = N3nExt(i+1,j+1,k ) + dnc%pq4(kk) * nit_dinc + N3nExt(i ,j ,k+1) = N3nExt(i ,j ,k+1) + dnc%pq5(kk) * nit_dinc + N3nExt(i+1,j ,k+1) = N3nExt(i+1,j ,k+1) + dnc%pq6(kk) * nit_dinc + N3nExt(i ,j+1,k+1) = N3nExt(i ,j+1,k+1) + dnc%pq7(kk) * nit_dinc + N3nExt(i+1,j+1,k+1) = N3nExt(i+1,j+1,k+1) + dnc%pq8(kk) * nit_dinc + + + endif + + enddo + +! we apply contribution in grd%variable + + slicevar(:,1:grd%km) = NutArray(1,:,:) + call ADD_PREVCORE_CONTRIB(N3nExt, grd%km, NutArray, slicevar) + ! call ADD_PREVCORE_CONTRIB(N3nExtended_3d, grd%km, grd%N3n_ad, grd%n3n_ad(1,:,:)) + + + +end subroutine v_densnut diff --git a/ver_hor_densnut.f90 b/ver_hor_densnut.f90 index be99ae2..3598586 100644 --- a/ver_hor_densnut.f90 +++ b/ver_hor_densnut.f90 @@ -1,4 +1,4 @@ -subroutine ver_hor_nut(NutArray, NutArrayAd, Var) +subroutine ver_hor_densnut(NutArray) !--------------------------------------------------------------------------- ! ! @@ -36,10 +36,11 @@ subroutine ver_hor_nut(NutArray, NutArrayAd, Var) use set_knd use grd_str - use eof_str +! use eof_str use cns_str use drv_str - use obs_str + use dnc_str +! use obs_str use mpi_str implicit none @@ -51,13 +52,13 @@ subroutine ver_hor_nut(NutArray, NutArrayAd, Var) type(DoubleGrid), allocatable, dimension(:) :: RecBuf1D(:) REAL(r8), allocatable, dimension(:,:,:) :: DefBufChl, DefBufChlAd REAL(r8) :: NutArray(grd%im,grd%jm,grd%km), NutArrayAd(grd%im,grd%jm,grd%km) - CHARACTER :: Var +! CHARACTER :: Var ione = 1 ! --- - ! Vertical EOFs - call veof_nut(NutArray, Var) + ! Apply density-nutrient increments and interp on vertical grid + call v_densnut(NutArray) !return ! goto 103 !No Vh @@ -248,4 +249,4 @@ subroutine ver_hor_nut(NutArray, NutArrayAd, Var) ! 103 continue -end subroutine ver_hor_nut +end subroutine ver_hor_densnut diff --git a/wrt_chl_stat.f90 b/wrt_chl_stat.f90 index 4c21db6..f2273eb 100644 --- a/wrt_chl_stat.f90 +++ b/wrt_chl_stat.f90 @@ -327,7 +327,7 @@ subroutine wrt_chl_stat - DEALLOCATE(DumpBio, ValuesToTest, MyConditions) + DEALLOCATE(DumpBio, ValuesToTest, MyConditions, LimitCorr) end subroutine wrt_chl_stat diff --git a/wrt_dia.f90 b/wrt_dia.f90 index e33276b..273da89 100644 --- a/wrt_dia.f90 +++ b/wrt_dia.f90 @@ -105,6 +105,18 @@ subroutine wrt_dia if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_att', status) endif + if(drv%dnc .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) + + 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 + status = nf90mpi_enddef(ncid) if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', status) @@ -169,26 +181,56 @@ subroutine wrt_dia 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 do j=1,grd%jm - do i=1,grd%im - if (drv%argo_obs .eq. 1) then - if (grd%msk(i,j,k) .eq. 1) then - DumpMatrix(i,j,k) = REAL(grd%o2o(i,j,k), 4) - else - DumpMatrix(i,j,k) = 1.e20 - endif + do i=1,grd%im + if (drv%argo_obs .eq. 1) then + if (grd%msk(i,j,k) .eq. 1) then + DumpMatrix(i,j,k) = REAL(grd%o2o(i,j,k), 4) else - DumpMatrix(i,j,k) = REAL(grd%o2o(i,j,k), 4 ) + DumpMatrix(i,j,k) = 1.e20 endif - enddo + else + DumpMatrix(i,j,k) = REAL(grd%o2o(i,j,k), 4 ) + endif + enddo enddo enddo status = nf90mpi_put_var_all(ncid,ido2o,DumpMatrix,MyStart,MyCount) if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all o2o', status) endif + + if(drv%dnc .eq. 1) then + do k=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) = REAL(grd%n3n(i,j,k), 4 ) + else + DumpMatrix(i,j,k) = 1.e20 + endif + enddo + enddo + enddo + status = nf90mpi_put_var_all(ncid,idn3n,DumpMatrix,MyStart,MyCount) + if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all n3n', status) + + do k=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) = REAL(grd%n3n(i,j,k)*bio%covn3n_n1p(i,j,k), 4 ) + else + DumpMatrix(i,j,k) = 1.e20 + endif + enddo + enddo + enddo + 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 status = nf90mpi_close(ncid) if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_close', status) diff --git a/wrt_nut_stat.f90 b/wrt_nut_stat.f90 index 91fed4a..f31ef86 100644 --- a/wrt_nut_stat.f90 +++ b/wrt_nut_stat.f90 @@ -66,7 +66,7 @@ subroutine wrt_nut_stat call MPI_Abort(Var3DCommunicator,-1,ierr) endif - if((bio%updateN1p .eq. 1) .and. (NNVar.lt.2)) then + if((bio%updateN1p .eq. 1) .and. (NNVar.lt.2) .and. (drv%dnc .eq. 0)) then write(*,*) "ERROR: Required phosphate update but NOT set in DA_params.f90" write(drv%dia,*) "ERROR: Required phosphate update but NOT set in DA_params.f90" call MPI_Barrier(Var3DCommunicator, ierr) @@ -97,7 +97,7 @@ subroutine wrt_nut_stat ! if the correction is negative, the correction must be reduced ! if(bio%n3n.eq.1) then ValuesToTest(i,j,k,1) = bio%InitialNut(i,j,k,1) + grd%n3n(i,j,k) - if(bio%updateN1p.eq.1) then + if((bio%updateN1p.eq.1) .or. (drv%dnc .eq. 1)) then ValuesToTest(i,j,k,2) = bio%InitialNut(i,j,k,2) + grd%n3n(i,j,k)*bio%covn3n_n1p(i,j,k) ! endif ! else From 423be85b9342b0950fb9bfa70dc6c2c5a70a90d5 Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Sun, 9 Nov 2025 22:41:57 +0100 Subject: [PATCH 05/13] Use of the correct parameter for densnut DA --- cnv_inn.f90 | 2 +- def_nml.f90 | 1 + def_nml_multi.f90 | 8 ++++---- drv_str.f90 | 3 +-- get_densincr_arg.f90 | 26 ++++++++++++++------------ get_obs.f90 | 2 +- int_par.f90 | 2 +- namelists/satfloat.20150101.nml | 3 +++ 8 files changed, 26 insertions(+), 21 deletions(-) diff --git a/cnv_inn.f90 b/cnv_inn.f90 index 83dcffc..e2df194 100644 --- a/cnv_inn.f90 +++ b/cnv_inn.f90 @@ -62,7 +62,7 @@ subroutine cnv_inn call ver_hor_nut(grd%n3n, grd%n3n_ad,'N') endif - if (drv%densnut .eq. 1) then + if (drv%dnc .eq. 1) then grd%n3n(:,:,:) = 0.0 call ver_hor_densnut(grd%n3n) endif diff --git a/def_nml.f90 b/def_nml.f90 index f99c065..7cf7efa 100644 --- a/def_nml.f90 +++ b/def_nml.f90 @@ -37,6 +37,7 @@ subroutine def_nml use ctl_str use mpi_str use bio_str + use dnc_str implicit none diff --git a/def_nml_multi.f90 b/def_nml_multi.f90 index 265a6cd..9cc06bb 100644 --- a/def_nml_multi.f90 +++ b/def_nml_multi.f90 @@ -47,14 +47,14 @@ subroutine def_nml_multi LOGICAL :: ApplyConditions INTEGER(i4) :: chl_assim, chl_upnut, nut, multiv, N3n, O2o, updateN1p - INTEGER(i4) :: nphyto, uniformL, anisL, densnut + INTEGER(i4) :: nphyto, uniformL, anisL, dnc 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, densnut, uniformL, anisL + NAMELIST /params/ sat_obs, argo, dnc, uniformL, anisL ! ------------------------------------------------------------------- @@ -119,7 +119,7 @@ subroutine def_nml_multi 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,*) ' Copuled density/nutrient ass densnut = ', densnut + write(drv%dia,*) ' Coupled density/nutrient ass dnc = ', dnc write(drv%dia,*) ' Set uniform correlation radius uniformL = ', uniformL write(drv%dia,*) ' Set anisotropy on corr radius anisL = ', anisL write(drv%dia,*) '------------------------------------------------------------' @@ -133,7 +133,7 @@ subroutine def_nml_multi drv%sat_obs = sat_obs drv%argo_obs = argo - drv%densnut = densnut + drv%dnc = dnc drv%uniformL = uniformL drv%anisL = anisL diff --git a/drv_str.f90 b/drv_str.f90 index 46d5492..f44a4a3 100644 --- a/drv_str.f90 +++ b/drv_str.f90 @@ -41,7 +41,7 @@ MODULE drv_str INTEGER(i4) :: MyCounter ! Number of iteration done by Tao solver 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) :: densnut ! Flag for the update of nutrients based on density increments + INTEGER(i4) :: dnc ! Flag for the update of nutrients based on density increments 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) @@ -49,7 +49,6 @@ MODULE drv_str INTEGER(i4) :: Verbose ! Flag for printing verbose output INTEGER(i4) :: nut ! Flag for nutrient assimilation INTEGER(i4) :: multiv ! Flag for multivariate assimilation - INTEGER(i4) :: dnc ! Flag for update of nitrate based on density increments END TYPE drv_t diff --git a/get_densincr_arg.f90 b/get_densincr_arg.f90 index 0ea307c..207605c 100644 --- a/get_densincr_arg.f90 +++ b/get_densincr_arg.f90 @@ -86,27 +86,29 @@ subroutine get_densincr_arg close (511) endif + call MPI_Bcast(TmpFlc, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) ! call MPI_Bcast(TmpPar, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpLon, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpLat, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpDpt, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) -! call MPI_Bcast(TmpTim, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) + ! call MPI_Bcast(TmpTim, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpInc, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpCorr, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpErr, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpStd, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - + ! Counting the number of observations that falls in the domain Counter = 0 do k=1,GlobalDncNum - 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(drv%dnc.eq.1) then - Counter = Counter + 1 - endif - endif - enddo + 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(drv%dnc.eq.1) then + Counter = Counter + 1 + endif + endif +enddo + if(drv%Verbose .eq. 1) & print*, "MyId", MyId, "has",Counter,"Density increments" @@ -374,10 +376,10 @@ subroutine int_par_dnc dnc%nc_global = 0 call MPI_Allreduce(dnc%nc, dnc%nc_global, 1, MPI_INT, MPI_SUM, Var3DCommunicator, ierr) - + if(MyId .eq. 0) then - write(drv%dia,*)'Real number of denisty increments: ',dnc%nc_global - print*,'Good denisty increments: ',dnc%nc_global + write(drv%dia,*)'Real number of density increments: ',dnc%nc_global + print*,'Good density increments: ',dnc%nc_global end if DEALLOCATE ( dnc%flg) diff --git a/get_obs.f90 b/get_obs.f90 index b506641..84b6687 100644 --- a/get_obs.f90 +++ b/get_obs.f90 @@ -52,7 +52,7 @@ subroutine get_obs call get_obs_sat ! Load density increments for Physical DA - if(drv%densnut .eq. 1) & + if(drv%dnc .eq. 1) & call get_densincr_arg end subroutine get_obs diff --git a/int_par.f90 b/int_par.f90 index 685d215..a5007fd 100644 --- a/int_par.f90 +++ b/int_par.f90 @@ -46,7 +46,7 @@ subroutine int_par call int_par_chl ! Load density increments for Physical DA - if(drv%densnut .eq. 1) & + if(drv%dnc .eq. 1) & call int_par_dnc end subroutine int_par diff --git a/namelists/satfloat.20150101.nml b/namelists/satfloat.20150101.nml index 5fbd536..12014b5 100644 --- a/namelists/satfloat.20150101.nml +++ b/namelists/satfloat.20150101.nml @@ -25,6 +25,7 @@ chl_assim = 1 chl_upnut = 0 nut = 1 + multiv = 0 nphyto = 4 chl_dep = 0.0 ncmp = 5 @@ -44,6 +45,7 @@ ApplyConditions = .true. ! 0-no satellite assimilation ! argo - 1-assimilate argo data ! - 0-no argo assimilation +! dnc - coupled density/nutrient assimilation ! uniformL - 1-non uniform radius ! - 0-uniform radius (rcf%L) ! anisL - 1-anisotropy @@ -53,6 +55,7 @@ ApplyConditions = .true. ¶ms sat_obs = 1 argo = 1 + dnc = 1 uniformL = 0 anisL = 0 / From 673875acd5b05257f6de2ab43b29089687251f8f Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Mon, 10 Nov 2025 14:46:22 +0100 Subject: [PATCH 06/13] Bug on allocation corrected, changes to better manage drv%dnc flag --- biovar.f90 | 8 ++++++-- def_cov.f90 | 8 +++----- sav_itr.f90 | 10 +++++----- wrt_nut_stat.f90 | 2 +- 4 files changed, 15 insertions(+), 13 deletions(-) diff --git a/biovar.f90 b/biovar.f90 index 6cf781e..1cc2271 100644 --- a/biovar.f90 +++ b/biovar.f90 @@ -115,7 +115,8 @@ subroutine biovar if (drv%chl_upnut .eq. 1) then call wrt_upd_nut else - call cp_nut_stat + if (drv%dnc .eq. 0) & + call cp_nut_stat endif endif endif @@ -130,7 +131,10 @@ subroutine biovar call wrt_nut_stat else if((bio%O2o .eq. 1) .and. (bio%N3n .eq. 0)) then call wrt_o2o_stat - call cp_nut_stat + if (drv%dnc .eq. 0) then + if(MyId .eq. 0) write(drv%dia,*) 'cp nut' + call cp_nut_stat + endif else if((bio%O2o .eq. 0) .and. (bio%N3n .eq. 1)) then call cp_o2o_stat call wrt_nut_stat diff --git a/def_cov.f90 b/def_cov.f90 index 0fd2d81..4939057 100644 --- a/def_cov.f90 +++ b/def_cov.f90 @@ -531,8 +531,10 @@ subroutine def_cov endif if(drv%nut.eq.1) then + if(MyId .eq. 0) write(drv%dia,*) 'read o2o and nut' call readNutStat - if(bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) then + if((bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) .or. (drv%dnc .eq. 1)) then + if(MyId .eq. 0) write(drv%dia,*) 'read nutcov' call readNutCov endif if(drv%chl_assim.eq.0) then @@ -547,8 +549,4 @@ subroutine def_cov call readNutCov endif - if(drv%dnc .eq. 1) then - call readNutStat - call readNutCov - endif end subroutine def_cov diff --git a/sav_itr.f90 b/sav_itr.f90 index f5c8e56..576c780 100644 --- a/sav_itr.f90 +++ b/sav_itr.f90 @@ -140,7 +140,7 @@ subroutine sav_itr endif if(drv%nut .eq. 1) then DEALLOCATE( bio%InitialNut) - if(bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) DEALLOCATE( bio%covn3n_n1p) + if((bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) .or. (drv%dnc .eq. 1)) DEALLOCATE(bio%covn3n_n1p) if(drv%chl_assim .eq. 0) then DEALLOCATE( bio%cquot, bio%pquot) DEALLOCATE( bio%InitialChl) !used in cp_chl_stat @@ -156,10 +156,10 @@ subroutine sav_itr if(bio%updateN1p.eq.1) DEALLOCATE( bio%covn3n_n1p) endif - if(drv%dnc .eq. 1) then - DEALLOCATE( bio%InitialNut) - DEALLOCATE( bio%covn3n_n1p) - endif + ! if(drv%dnc .eq. 1) then + ! DEALLOCATE( bio%InitialNut) + ! DEALLOCATE( bio%covn3n_n1p) + ! endif DEALLOCATE(SurfaceWaterPoints) diff --git a/wrt_nut_stat.f90 b/wrt_nut_stat.f90 index f31ef86..cf82cc2 100644 --- a/wrt_nut_stat.f90 +++ b/wrt_nut_stat.f90 @@ -59,7 +59,7 @@ subroutine wrt_nut_stat TimeArr(1) = DA_JulianDate - if(bio%n3n .eq. 0) then + if((bio%n3n .eq. 0) .and. (drv%dnc .eq. 0)) then write(*,*) "ERROR: Nitrate to be assimilated NOT set in namelist" write(drv%dia,*) "ERROR: Nitrate to be assimilated NOT set in namelist" call MPI_Barrier(Var3DCommunicator, ierr) From 2754f2cb3dfc8f895a1ce27a93685e5c881d084f Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Mon, 9 Feb 2026 15:59:17 +0100 Subject: [PATCH 07/13] FIRST!!! Version with density-nutrient EOFs To be compiled and tested! --- biovar.f90 | 18 ++---- clean_mem.f90 | 8 +-- cnv_inn.f90 | 9 ++- costf.f90 | 10 ++- def_cov.f90 | 10 ++- def_nml.f90 | 2 + def_nml_multi.f90 | 4 +- dnc_str.f90 | 6 +- eof_str.f90 | 4 ++ filename_mod.f90 | 4 +- get_densincr_arg.f90 | 40 ++++++------ mpi_str.f90 | 2 +- obs_arg.f90 | 2 +- obs_dnc.f90 | 81 ++++++++++++++++++++++++ obs_dnc_ad.f90 | 96 +++++++++++++++++++++++++++++ obs_vec.f90 | 16 ++++- obsop.f90 | 4 ++ obsop_ad.f90 | 3 + rdeofs_dnc.f90 | 144 +++++++++++++++++++++++++++++++++++++++++++ readGrid.f90 | 14 +++-- res_inc.f90 | 4 +- resid.f90 | 13 ++++ veof_dnc.f90 | 81 ++++++++++++++++++++++++ veof_dnc_ad.f90 | 128 ++++++++++++++++++++++++++++++++++++++ veof_nut.f90 | 30 +++++---- ver_hor_nut.f90 | 6 +- ver_hor_nut_ad.f90 | 41 ++++++------ 27 files changed, 690 insertions(+), 90 deletions(-) create mode 100644 obs_dnc.f90 create mode 100644 obs_dnc_ad.f90 create mode 100644 rdeofs_dnc.f90 create mode 100644 veof_dnc.f90 create mode 100644 veof_dnc_ad.f90 diff --git a/biovar.f90 b/biovar.f90 index 1cc2271..667b67f 100644 --- a/biovar.f90 +++ b/biovar.f90 @@ -112,11 +112,9 @@ subroutine biovar ! 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 call cp_o2o_stat + call cp_nut_stat if (drv%chl_upnut .eq. 1) then call wrt_upd_nut - else - if (drv%dnc .eq. 0) & - call cp_nut_stat endif endif endif @@ -126,24 +124,18 @@ subroutine biovar endif if ((drv%nut .eq. 1) .or. (drv%multiv.eq.1)) then - if((bio%O2o .eq. 1) .and. (bio%N3n .eq. 1)) then + if((bio%O2o .eq. 1) .and. ((bio%N3n .eq. 1) .or. (drv%dnc .eq. 1))) then call wrt_o2o_stat call wrt_nut_stat - else if((bio%O2o .eq. 1) .and. (bio%N3n .eq. 0)) then + else if((bio%O2o .eq. 1) .and. (bio%N3n .eq. 0) .and. (drv%dnc .eq. 0)) then call wrt_o2o_stat - if (drv%dnc .eq. 0) then - if(MyId .eq. 0) write(drv%dia,*) 'cp nut' - call cp_nut_stat - endif - else if((bio%O2o .eq. 0) .and. (bio%N3n .eq. 1)) then + call cp_nut_stat + else if((bio%O2o .eq. 0) .and. ((bio%N3n .eq. 1) .or. (drv%dnc .eq. 1))) then call cp_o2o_stat call wrt_nut_stat endif endif - if (drv%dnc .eq. 1) then - call wrt_nut_stat - endif call sav_itr if(MyId .eq. 0) write(drv%dia,*) 'out of sav_itr ' diff --git a/clean_mem.f90 b/clean_mem.f90 index 355735d..e1a0d56 100644 --- a/clean_mem.f90 +++ b/clean_mem.f90 @@ -72,10 +72,10 @@ subroutine clean_mem ! density increments if(drv%dnc .eq. 1) then DEALLOCATE ( dnc%flc) - DEALLOCATE ( dnc%inc) - DEALLOCATE ( dnc%corr) + DEALLOCATE ( dnc%res) + ! DEALLOCATE ( dnc%corr) DEALLOCATE ( dnc%err) - DEALLOCATE ( dnc%std) + ! DEALLOCATE ( dnc%std) DEALLOCATE ( dnc%ib, dnc%jb, dnc%kb) DEALLOCATE ( dnc%pq1, dnc%pq2, dnc%pq3, dnc%pq4) DEALLOCATE ( dnc%pq5, dnc%pq6, dnc%pq7, dnc%pq8) @@ -91,7 +91,7 @@ subroutine clean_mem DEALLOCATE(RecCountX3D_chl, RecDisplX3D_chl) DEALLOCATE(ChlExtended) - DEALLOCATE(ChlExtended_3d,N3nExtended_3d,O2oExtended_3d) + DEALLOCATE(ChlExtended_3d,N3nExtended_3d,O2oExtended_3d,DncExtended_3d) DEALLOCATE(SendBottom, RecTop) DEALLOCATE(SendTop, RecBottom) DEALLOCATE(SendTop_2d, RecBottom_2d) diff --git a/cnv_inn.f90 b/cnv_inn.f90 index e2df194..6eb4fdf 100644 --- a/cnv_inn.f90 +++ b/cnv_inn.f90 @@ -48,8 +48,11 @@ subroutine cnv_inn call ver_hor_chl endif if(drv%nut .eq. 1) then - if(bio%N3n .eq. 1) then + if((bio%N3n .eq. 1) .or. (drv%dnc .eq. 1)) then call ver_hor_nut(grd%n3n, grd%n3n_ad,'N') + ! if (drv%dnc .eq. 1) then + ! call ver_hor_nut(grd%n3n) + ! endif endif if(bio%O2o .eq. 1) then call ver_hor_nut(grd%o2o, grd%o2o_ad,'O') @@ -62,10 +65,6 @@ subroutine cnv_inn call ver_hor_nut(grd%n3n, grd%n3n_ad,'N') endif - if (drv%dnc .eq. 1) then - grd%n3n(:,:,:) = 0.0 - call ver_hor_densnut(grd%n3n) - endif ! --- ! Apply biological repartition of the chlorophyll if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) & diff --git a/costf.f90 b/costf.f90 index 41570c1..3f69a1e 100644 --- a/costf.f90 +++ b/costf.f90 @@ -61,8 +61,11 @@ subroutine costf call ver_hor_chl endif if(drv%nut .eq. 1) then - if(bio%N3n .eq. 1) then + if((bio%N3n .eq. 1) .or. (drv%dnc .eq. 1)) then call ver_hor_nut(grd%n3n, grd%n3n_ad, 'N') + if(drv%dnc .eq. 1) then + call ver_hor_nut(grd%dnc,grd%dnc_ad,'D') + endif endif if(bio%O2o .eq. 1) then call ver_hor_nut(grd%o2o, grd%o2o_ad, 'O') @@ -125,8 +128,11 @@ subroutine costf call ver_hor_chl_ad endif if(drv%nut .eq. 1) then - if(bio%N3n .eq. 1) then + if((bio%N3n .eq. 1) .or. (drv%dnc .eq. 1)) then call ver_hor_nut_ad(grd%n3n, grd%n3n_ad, 'N') + if(drv%dnc .eq. 1) then + call ver_hor_nut_ad(grd%dnc, grd%dnc_ad, 'D') + endif endif if(bio%O2o .eq. 1) then call ver_hor_nut_ad(grd%o2o, grd%o2o_ad, 'O') diff --git a/def_cov.f90 b/def_cov.f90 index 4939057..7217c62 100644 --- a/def_cov.f90 +++ b/def_cov.f90 @@ -465,12 +465,17 @@ subroutine def_cov else ros%neof_n3n = 0 endif + if(drv%dnc .eq. 1) then + call rdeofs_dnc + else + ros%neof_dnc = 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 + ros%neof_nut = ros%neof_n3n + ros%neof_o2o + ros%neof_dnc else ros%neof_nut = 0 endif @@ -479,6 +484,7 @@ subroutine def_cov ros%neof_chl = 0 ros%neof_n3n = 0 ros%neof_o2o = 0 + ros%neof_dnc = 0 call rdeofs_multi endif @@ -533,7 +539,7 @@ subroutine def_cov if(drv%nut.eq.1) then if(MyId .eq. 0) write(drv%dia,*) 'read o2o and nut' call readNutStat - if((bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) .or. (drv%dnc .eq. 1)) then + if((bio%N3n.eq.1 .or. drv%dnc .eq. 1) .AND. bio%updateN1p.eq.1) then if(MyId .eq. 0) write(drv%dia,*) 'read nutcov' call readNutCov endif diff --git a/def_nml.f90 b/def_nml.f90 index 7cf7efa..4c92e1e 100644 --- a/def_nml.f90 +++ b/def_nml.f90 @@ -100,6 +100,7 @@ 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 dens-nut EOFs: neof_dnc = ', neof_dnc 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 @@ -116,6 +117,7 @@ subroutine def_nml ros%neof_chl = neof_chl ros%neof_n3n = neof_n3n ros%neof_o2o = neof_o2o + ros%neof_dnc = neof_dnc ros%neof_multi = neof_multi ros%kmchl = kmchl ros%kmnit = kmnit diff --git a/def_nml_multi.f90 b/def_nml_multi.f90 index 9cc06bb..034be84 100644 --- a/def_nml_multi.f90 +++ b/def_nml_multi.f90 @@ -84,7 +84,7 @@ subroutine def_nml_multi 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,*) ' Nutrient assimilation ALSO for dnc 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 @@ -119,7 +119,7 @@ subroutine def_nml_multi 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,*) ' Coupled density/nutrient ass dnc = ', dnc + write(drv%dia,*) ' Coupled density/nutrient DA dnc = ', dnc write(drv%dia,*) ' Set uniform correlation radius uniformL = ', uniformL write(drv%dia,*) ' Set anisotropy on corr radius anisL = ', anisL write(drv%dia,*) '------------------------------------------------------------' diff --git a/dnc_str.f90 b/dnc_str.f90 index da66cd6..ead588c 100644 --- a/dnc_str.f90 +++ b/dnc_str.f90 @@ -52,10 +52,10 @@ MODULE dnc_str REAL(r8), POINTER :: dpt(:) ! Depth ! REAL(r8), POINTER :: tim(:) ! Time REAL(r8), POINTER :: inc(:) ! Increments - REAL(r8), POINTER :: corr(:) ! Correlations + ! REAL(r8), POINTER :: corr(:) ! Correlations REAL(r8), POINTER :: err(:) ! Nitrate std (error) - REAL(r8), POINTER :: std(:) ! Density std - ! REAL(r8), POINTER :: res(:) ! residual + ! REAL(r8), POINTER :: std(:) ! Density std + 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 INTEGER(i8), POINTER :: jb(:) ! j index of the nearest south point diff --git a/eof_str.f90 b/eof_str.f90 index 75cd61c..0d1c82e 100644 --- a/eof_str.f90 +++ b/eof_str.f90 @@ -43,6 +43,7 @@ MODULE eof_str 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) :: neof_dnc ! No. of EOFs for coupled dens-nit 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 @@ -62,6 +63,9 @@ MODULE eof_str REAL(r8), POINTER :: eva_o2o(:,:) ! Eigenvalues REAL(r8), POINTER :: evc_multi(:,:,:) ! Eigenvectors REAL(r8), POINTER :: eva_multi(:,:) ! Eigenvalues + ! Optional EOFs for density-nutrient covariance + REAL(r8), POINTER :: evc_dnc(:,:,:) ! Eigenvectors for density-nutrient covariance + REAL(r8), POINTER :: eva_dnc(:,:) ! Eigenvalues for density-nutrient covariance #endif diff --git a/filename_mod.f90 b/filename_mod.f90 index a91a9ce..cf6bcb6 100644 --- a/filename_mod.f90 +++ b/filename_mod.f90 @@ -6,6 +6,7 @@ MODULE FILENAMES 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) :: EOF_FILE_DNC != 'eofs_dnc.nc' character (LEN=1024) :: STD_FILE_MULTI != 'std_multi.nc' character (LEN=1024) :: MISFIT_FILE != 'chl_mis.nc' character (LEN=1024) :: NUTCOV_FILE != 'crosscorrs.nc' @@ -30,6 +31,7 @@ SUBROUTINE SETFILENAMES EOF_FILE_N3N = 'eofs_n3n.nc' EOF_FILE_O2O = 'eofs_o2o.nc' EOF_FILE_MULTI = 'eofs_multi.nc' +EOF_FILE_DNC = 'eofs_dnc.nc' STD_FILE_MULTI = 'std_multi.nc' MISFIT_FILE = 'chl_mis.nc' NUTCOV_FILE = 'crosscorrs.nc' @@ -41,7 +43,7 @@ SUBROUTINE SETFILENAMES RCORR_FILE = 'chl_rad_corr.nc' ARGO_FILE = 'arg_mis.dat' ANIS_FILE = 'gradsal.nc' - INCR_FILE = 'density_increments.dat' + INCR_FILE = 'density_increments.nc' END SUBROUTINE SETFILENAMES diff --git a/get_densincr_arg.f90 b/get_densincr_arg.f90 index 207605c..6d02b2a 100644 --- a/get_densincr_arg.f90 +++ b/get_densincr_arg.f90 @@ -42,8 +42,8 @@ subroutine get_densincr_arg INTEGER(i4) :: k INTEGER(i4) :: i1, kk, i REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpFlc, TmpLon, TmpLat - REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpDpt, TmpErr, TmpStd - REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpInc, TmpCorr + REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpDpt, TmpErr!, TmpStd + REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpInc!, TmpCorr INTEGER(i4) :: GlobalDncNum, Counter, ierr character(len=1024) :: filename @@ -70,8 +70,8 @@ subroutine get_densincr_arg ALLOCATE( TmpFlc(GlobalDncNum)) !, TmpPar(GlobalDncNum)) ALLOCATE( TmpLon(GlobalDncNum), TmpLat(GlobalDncNum)) ALLOCATE( TmpDpt(GlobalDncNum))!, TmpTim(GlobalDncNum)) - ALLOCATE( TmpInc(GlobalDncNum), TmpCorr(GlobalDncNum)) - ALLOCATE( TmpErr(GlobalDncNum), TmpStd(GlobalDncNum)) + ALLOCATE( TmpInc(GlobalDncNum))!, TmpCorr(GlobalDncNum)) + ALLOCATE( TmpErr(GlobalDncNum))!, TmpStd(GlobalDncNum)) if(MyId .eq. 0) then ! process 0 reads all the density increments @@ -80,8 +80,8 @@ subroutine get_densincr_arg TmpFlc(k), & !TmpPar(k), & TmpLon(k), TmpLat(k), & TmpDpt(k), &!TmpTim(k), & - TmpInc(k), TmpCorr(k), & - TmpErr(k), TmpStd(k) + TmpInc(k), &!TmpCorr(k), & + TmpErr(k) !, TmpStd(k) end do close (511) endif @@ -94,18 +94,18 @@ subroutine get_densincr_arg call MPI_Bcast(TmpDpt, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) ! call MPI_Bcast(TmpTim, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpInc, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpCorr, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) +! call MPI_Bcast(TmpCorr, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) call MPI_Bcast(TmpErr, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) - call MPI_Bcast(TmpStd, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) +! call MPI_Bcast(TmpStd, GlobalDncNum, MPI_REAL8, 0, Var3DCommunicator, ierr) ! Counting the number of observations that falls in the domain Counter = 0 do k=1,GlobalDncNum 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(drv%dnc.eq.1) then + ! if(drv%dnc.eq.1) then Counter = Counter + 1 - endif + ! endif endif enddo @@ -118,9 +118,10 @@ subroutine get_densincr_arg ALLOCATE ( dnc%flg(dnc%no), dnc%flc(dnc%no)) !, dnc%par(dnc%no)) ALLOCATE ( dnc%lon(dnc%no), dnc%lat(dnc%no), dnc%dpt(dnc%no)) !, dnc%tim(dnc%no)) ALLOCATE ( dnc%inc(dnc%no)) - ALLOCATE ( dnc%corr(dnc%no)) +! ALLOCATE ( dnc%corr(dnc%no)) ALLOCATE ( dnc%err(dnc%no)) - ALLOCATE ( dnc%std(dnc%no)) + ALLOCATE ( dnc%res(dnc%no)) +! ALLOCATE ( dnc%std(dnc%no)) ALLOCATE ( dnc%ib(dnc%no), dnc%jb(dnc%no), dnc%kb(dnc%no)) ALLOCATE ( dnc%pb(dnc%no), dnc%qb(dnc%no), dnc%rb(dnc%no)) ALLOCATE ( dnc%pq1(dnc%no), dnc%pq2(dnc%no), dnc%pq3(dnc%no), dnc%pq4(dnc%no)) @@ -130,19 +131,19 @@ subroutine get_densincr_arg do k=1,GlobalDncNum 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(drv%dnc.eq.1) then + ! if(drv%dnc.eq.1) then Counter = Counter + 1 dnc%flc(Counter) = TmpFlc(k) ! dnc%par(Counter) = TmpPar(k) dnc%lon(Counter) = TmpLon(k) dnc%lat(Counter) = TmpLat(k) dnc%dpt(Counter) = TmpDpt(k) - dnc%inc(Counter) = TmpInc(k) - dnc%corr(Counter) = TmpCorr(k) + dnc%res(Counter) = TmpInc(k) !called res for analogy with get_obs_arg + ! dnc%corr(Counter) = TmpCorr(k) dnc%err(Counter) = TmpErr(k) - dnc%std(Counter) = TmpStd(k) + ! dnc%std(Counter) = TmpStd(k) ! dnc%ino(Counter) = TmpIno(k) - endif + ! endif endif enddo @@ -177,6 +178,7 @@ subroutine get_densincr_arg if(dnc%flg(k).eq.1)then dnc%nc = dnc%nc + 1 else + dnc%res(k) = 0. dnc%inc(k) = 0. dnc%pq1(k) = 0. dnc%pq2(k) = 0. @@ -194,9 +196,9 @@ subroutine get_densincr_arg DEALLOCATE( TmpLon, TmpLat) DEALLOCATE( TmpDpt)!, TmpTim) DEALLOCATE( TmpInc)!, TmpTim) - DEALLOCATE( TmpCorr)!, TmpTim) +! DEALLOCATE( TmpCorr)!, TmpTim) DEALLOCATE( TmpErr) - DEALLOCATE( TMPStd) +! DEALLOCATE( TMPStd) ! DEALLOCATE( TmpIno) end subroutine get_densincr_arg diff --git a/mpi_str.f90 b/mpi_str.f90 index 7858202..4816d24 100644 --- a/mpi_str.f90 +++ b/mpi_str.f90 @@ -62,7 +62,7 @@ MODULE mpi_str 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 + REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: ChlExtended_3d, N3nExtended_3d, O2oExtended_3d, DncExtended_3d CONTAINS diff --git a/obs_arg.f90 b/obs_arg.f90 index 39d95f8..ace2e80 100644 --- a/obs_arg.f90 +++ b/obs_arg.f90 @@ -52,7 +52,7 @@ subroutine obs_arg 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 + 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 diff --git a/obs_dnc.f90 b/obs_dnc.f90 new file mode 100644 index 0000000..4d6c5cd --- /dev/null +++ b/obs_dnc.f90 @@ -0,0 +1,81 @@ +subroutine obs_dnc + + + !--------------------------------------------------------------------------- + ! ! + ! 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 . ! + ! ! + !--------------------------------------------------------------------------- + + !----------------------------------------------------------------------- + ! ! + ! Apply observational operator for ARGO floats ! + ! ! + ! Version 1: S.Dobricic 2006 ! + !----------------------------------------------------------------------- + + + use set_knd + use grd_str + use eof_str + use dnc_str + use mpi_str + use drv_str + use bio_str + + implicit none + + INTEGER(i4) :: i, j, k, kk, condc, condn, my_km + + 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. drv%dnc.eq.1 ) .or. (drv%multiv.eq.1)) then + ! condn = 1 + call EXTEND_2D( grd%dnc, grd%km, DncExtended_3d ) + ! endif + ! if (bio%O2o.eq.1 ) & + ! call EXTEND_2D( grd%O2o, grd%km, O2oExtended_3d ) + + + + do kk = 1,dnc%no + if(dnc%flc(kk).eq.1) then + + + dnc%inc(kk) = & + dnc%pq1(kk) * DncExtended_3d(i ,j ,k) + & + dnc%pq2(kk) * DncExtended_3d(i+1,j ,k ) + & + dnc%pq3(kk) * DncExtended_3d(i ,j+1,k ) + & + dnc%pq4(kk) * DncExtended_3d(i+1,j+1,k ) + & + dnc%pq5(kk) * DncExtended_3d(i ,j ,k+1) + & + dnc%pq6(kk) * DncExtended_3d(i+1,j ,k+1) + & + dnc%pq7(kk) * DncExtended_3d(i ,j+1,k+1) + & + dnc%pq8(kk) * DncExtended_3d(i+1,j+1,k+1) + endif + enddo + + +end subroutine obs_dnc diff --git a/obs_dnc_ad.f90 b/obs_dnc_ad.f90 new file mode 100644 index 0000000..c55a01b --- /dev/null +++ b/obs_dnc_ad.f90 @@ -0,0 +1,96 @@ +subroutine obs_dnc_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 . ! + ! ! + !--------------------------------------------------------------------------- + + !----------------------------------------------------------------------- + ! ! + ! Apply observational operator for density (adjoint) ! + ! ! + ! Version 1: S.Dobricic 2006 ! + !----------------------------------------------------------------------- + + + use set_knd + use grd_str + use eof_str + use obs_str + use mpi_str + use filenames + use drv_str + use bio_str + + implicit none + + 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 ) + call EXTEND_2D( grd%dnc_ad, grd%km, DncExtended_3d ) + + + do kk = 1,dnc%no + + i=dnc%ib(kk) + j=dnc%jb(kk) + k=dnc%kb(kk) + + if(dnc%flc(kk).eq.1)then + + obs%k = obs%k + 1 + obsg = obs%gra(obs%k) + + DncExtended_3d(i ,j ,k ) = DncExtended_3d(i ,j ,k ) + dnc%pq1(kk) * obsg + DncExtended_3d(i+1,j ,k ) = DncExtended_3d(i+1,j ,k ) + dnc%pq2(kk) * obsg + DncExtended_3d(i ,j+1,k ) = DncExtended_3d(i ,j+1,k ) + dnc%pq3(kk) * obsg + DncExtended_3d(i+1,j+1,k ) = DncExtended_3d(i+1,j+1,k ) + dnc%pq4(kk) * obsg + DncExtended_3d(i ,j ,k+1) = DncExtended_3d(i ,j ,k+1) + dnc%pq5(kk) * obsg + DncExtended_3d(i+1,j ,k+1) = DncExtended_3d(i+1,j ,k+1) + dnc%pq6(kk) * obsg + DncExtended_3d(i ,j+1,k+1) = DncExtended_3d(i ,j+1,k+1) + dnc%pq7(kk) * obsg + DncExtended_3d(i+1,j+1,k+1) = DncExtended_3d(i+1,j+1,k+1) + dnc%pq8(kk) * obsg + endif + + enddo + + + slicevar(:,1:my_km) = grd%dnc_ad(1,:,1:my_km) + call ADD_PREVCORE_CONTRIB(DncExtended_3d, my_km, grd%dnc_ad, slicevar(:,1:my_km)) + + + +end subroutine obs_dnc_ad diff --git a/obs_vec.f90 b/obs_vec.f90 index dd9d7fb..d8e6363 100644 --- a/obs_vec.f90 +++ b/obs_vec.f90 @@ -34,6 +34,7 @@ subroutine obs_vec use drv_str use obs_str use mpi_str + use dnc_str implicit none @@ -42,7 +43,7 @@ subroutine obs_vec ! ------- ! Define observational vector - obs%no = sat%nc + arg%nc + obs%no = sat%nc + arg%nc + dnc%nc if(MyId .eq. 0) & write(drv%dia,*) ' Total number of observations: ', obs%no @@ -84,5 +85,18 @@ subroutine obs_vec endif + ! Observations of density increments + if(drv%dnc .eq. 1) then + do i=1,dnc%no + if(dnc%flc(i).eq.1)then + k=k+1 + obs%res(k) = dnc%res(i) + obs%err(k) = dnc%err(i) + endif + enddo + + DEALLOCATE(dnc%res, dnc%err) + + endif end subroutine obs_vec diff --git a/obsop.f90 b/obsop.f90 index 466e1f1..180ee4b 100644 --- a/obsop.f90 +++ b/obsop.f90 @@ -55,6 +55,10 @@ subroutine obsop if(drv%sat_obs .eq. 1) & call obs_sat + ! Density increments + if(drv%dnc .eq. 1) & + call obs_dnc + call MPI_Barrier(Var3DCommunicator, ierr) end subroutine obsop diff --git a/obsop_ad.f90 b/obsop_ad.f90 index cc90231..f455d9c 100644 --- a/obsop_ad.f90 +++ b/obsop_ad.f90 @@ -50,6 +50,9 @@ subroutine obsop_ad if(drv%sat_obs .eq. 1) & call obs_sat_ad + ! Density increments + if(drv%dnc .eq. 1) & + call obs_dnc_ad ! --- ! Apply biological repartition of the chlorophyll if((drv%chl_assim .eq. 1) .or. (drv%multiv .eq. 1)) & diff --git a/rdeofs_dnc.f90 b/rdeofs_dnc.f90 new file mode 100644 index 0000000..3be64a6 --- /dev/null +++ b/rdeofs_dnc.f90 @@ -0,0 +1,144 @@ +subroutine rdeofs_dnc + + !--------------------------------------------------------------------------- + ! ! + ! 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) :: neofs, nlevs, nregs + integer(KIND=MPI_OFFSET_KIND) :: GlobalStart(3), GlobalCount(3) + real(4), allocatable :: x3(:,:,:), x2(:,:) + + + stat = nf90mpi_open(Var3DCommunicator, trim(EOF_FILE_DNC), NF90_NOWRITE, MPI_INFO_NULL, ncid) + if (stat /= nf90_noerr) call handle_err("nf90mpi_open "//trim(EOF_FILE_DNC), 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 nut-density are: ',ros%nreg, 2*ros%kmt, neofs + write(drv%dia,*)'Uses ',ros%neof_dnc,' 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_dnc .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_dnc .lt. neofs) then + + if(MyId .eq. 0) then + write(drv%dia,*)'Warning: ros%neof_dnc < neofs!' + write(drv%dia,*)'ros%neof_dnc =', ros%neof_dnc + write(drv%dia,*)'neofs =', neofs + write(drv%dia,*)'continue using ros%neof_dnc' + write(*,*)'Warning: ros%neof_dnc < neofs!' + write(*,*)'ros%neof_dnc =', ros%neof_dnc + write(*,*)'neofs =', neofs + write(*,*)'continue using ros%neof_dnc' + endif + endif + + if(2*ros%kmt .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_dnc( ros%nreg, 2*ros%kmt, ros%neof_dnc) ) ; ros%evc_dnc = huge(ros%evc_dnc(1,1,1)) + ALLOCATE ( ros%eva_dnc( ros%nreg, ros%neof_dnc) ) ; ros%eva_dnc = huge(ros%eva_dnc(1,1)) + ALLOCATE ( x3( ros%nreg, 2*ros%kmt, ros%neof_dnc) ) + ALLOCATE ( x2( ros%nreg, ros%neof_dnc) ) + GlobalStart(:) = 1 + GlobalCount(1) = ros%nreg + GlobalCount(2) = 2*ros%kmt + GlobalCount(3) = ros%neof_dnc + + 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_dnc(:,:,:) = x3(:,:,:) + + GlobalCount(1) = ros%nreg + GlobalCount(2) = ros%neof_dnc + + 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_dnc(:,:) = x2(:,:) + + ! DECOMMENT FOLLOWING TWO LINES TO MAKE FILTER TEST + ! ros%evc(:,:,:) = 1. + ! ros%eva(:,:) = 1. + + stat = nf90mpi_close(ncid) + if (stat /= nf90_noerr) call handle_err("nf90mpi_close", stat) + + DEALLOCATE(x3, x2) + +end subroutine rdeofs_dnc + + diff --git a/readGrid.f90 b/readGrid.f90 index 25ff1f3..2290f4c 100644 --- a/readGrid.f90 +++ b/readGrid.f90 @@ -111,6 +111,7 @@ subroutine readGrid end if ALLOCATE(N3nExtended_3d (grd%im+1, grd%jm, grd%km)) ALLOCATE(O2oExtended_3d (grd%im+1, grd%jm, grd%km)) + ALLOCATE(DncExtended_3d (grd%im+1, grd%jm, grd%km)) @@ -151,6 +152,12 @@ subroutine readGrid 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 + if(drv%dnc .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)) + ALLOCATE ( grd%dnc(grd%im,grd%jm,grd%km) ) ; grd%dnc = huge(grd%dnc(1,1,1)) + ALLOCATE ( grd%dnc_ad(grd%im,grd%jm,grd%km) ) ; grd%dnc_ad = huge(grd%dnc_ad(1,1,1)) + endif endif else if(drv%multiv .eq. 1) then @@ -162,15 +169,12 @@ subroutine readGrid 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 - if (drv%dnc .eq. 1) then - ALLOCATE (grd%n3n(grd%im,grd%jm,grd%km) ) ; grd%n3n = huge(grd%n3n(1,1,1)) - endif 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)) ALLOCATE ( x1(grd%km) ) ; x1 = huge(x1(1)) - if (drv%argo_obs .eq. 1) then + if ((drv%argo_obs .eq. 1) .or. (drv%dnc)) then ALLOCATE ( grd%lon(grd%im,grd%jm)) ; grd%lon = huge(grd%lon(1,1)) ALLOCATE ( grd%lat(grd%im,grd%jm)) ; grd%lat = huge(grd%lat(1,1)) endif @@ -187,7 +191,7 @@ subroutine readGrid if (ierr .ne. NF90_NOERR ) call handle_err('nfmpi_get_vara_real_all dy', ierr) grd%dy(:,:) = x2(:,:) - if (drv%argo_obs .eq. 1) then + if ((drv%argo_obs .eq. 1) .or. (drv%dnc)) then ierr = nf90mpi_inq_varid (ncid, 'lon', VarId) if (ierr .ne. NF90_NOERR ) call handle_err('nf90mpi_inq_varid', ierr) ierr = nfmpi_get_vara_real_all (ncid, VarId, MyStart, MyCount, x2) diff --git a/res_inc.f90 b/res_inc.f90 index edc891b..7c99611 100644 --- a/res_inc.f90 +++ b/res_inc.f90 @@ -43,8 +43,10 @@ subroutine res_inc end if if (drv%nut .eq. 1) then - if (bio%n3n .eq. 1) & + if ((bio%n3n .eq. 1) .or. (drv%dnc .eq. 1)) & grd%n3n_ad(:,:,:) = 0.0 + if (drv%dnc .eq. 1) & + grd%dnc(:,:,:) = 0.0 if (bio%o2o .eq. 1) & grd%o2o_ad(:,:,:) = 0.0 endif diff --git a/resid.f90 b/resid.f90 index faba056..3253024 100644 --- a/resid.f90 +++ b/resid.f90 @@ -64,4 +64,17 @@ subroutine resid enddo endif + + ! --- + ! Density increments + if(drv%dnc .eq. 1) then + do i=1,dnc%no + if(dnc%flc(i).eq.1)then + k = k + 1 + obs%inc(k) = dnc%inc(i) + obs%amo(k) = ( obs%inc(k) - obs%res(k) ) / obs%err(k) + endif + enddo + endif + end subroutine resid diff --git a/veof_dnc.f90 b/veof_dnc.f90 new file mode 100644 index 0000000..26aa304 --- /dev/null +++ b/veof_dnc.f90 @@ -0,0 +1,81 @@ +subroutine veof_dnc +!anna +!--------------------------------------------------------------------------- +! ! +! 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 +! ! +! 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, MyNEofs, ierr + REAL(r8), DIMENSION ( grd%im, grd%jm) :: egm + REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: eva + REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: evc + + my_km = grd%km + MyNEofs = ros%neof_dnc + offset = ros%neof_chl + ros%neof_n3n + + + ALLOCATE (eva(ros%nreg,MyNEofs)); eva = huge(eva(1,1)) + ALLOCATE (evc(ros%nreg,my_km,MyNEofs)); evc = huge(evc(1,1,1)) + + eva(:,:) = ros%eva_dnc(:,:) + evc(:,1:my_km,:) = ros%evc_dnc(:,my_km+1:my_km*2,:) + + grd%dnc(:,:,:) = 0.0 + + !cdir noconcur + do n=1,MyNEofs + + egm(:,:) = 0.0 + + do j=1,grd%jm + do i=1,grd%im + egm(i,j) = eva(grd%reg(i,j),n) * grd%ro( i, j, n+offset) + enddo + enddo + + ! 3D variables + do k=1,my_km ! OMP + do j=1,grd%jm + do i=1,grd%im + grd%dnc(i,j,k) = grd%dnc(i,j,k) + evc(grd%reg(i,j),k,n) * egm(i,j) + enddo + enddo + enddo + enddo + + DEALLOCATE(eva,evc) + +end subroutine veof_dnc diff --git a/veof_dnc_ad.f90 b/veof_dnc_ad.f90 new file mode 100644 index 0000000..6f57767 --- /dev/null +++ b/veof_dnc_ad.f90 @@ -0,0 +1,128 @@ +subroutine veof_dnc_ad(NutArrayAd, Var) + +!--------------------------------------------------------------------------- +! ! +! 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 + + implicit none + + 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 + MyNEofs = ros%neof_dnc + 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)) + + eva = ros%eva_dnc + evc = ros%evc_dnc + ! if((drv%nut .eq.1) .and. (drv%multiv .eq. 0)) then + ! if(Var .eq. 'N') then + ! 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 + grd%ro_ad(:,:,n+offset) = 0.0 ! OMP + enddo + +!$OMP PARALLEL & +!$OMP PRIVATE(i,j,k,k1,n) & +!$OMP PRIVATE(egm) +!$OMP DO + do n=1,MyNEofs + + egm(:,:) = 0.0 + + ! 3D variables + + do k=1,my_km ! OMP + do j=1,grd%jm + do i=1,grd%im + egm(i,j) = egm(i,j) + evc(grd%reg(i,j), k,n) * NutArrayAd(i,j,k) + enddo + enddo + enddo + + + do j=1,grd%jm + do i=1,grd%im + 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) + enddo + enddo + +enddo +!$OMP END DO +!$OMP END PARALLEL + +DEALLOCATE(eva,evc) + +end subroutine veof_dnc_ad diff --git a/veof_nut.f90 b/veof_nut.f90 index 49fe217..abff91b 100644 --- a/veof_nut.f90 +++ b/veof_nut.f90 @@ -42,7 +42,7 @@ subroutine veof_nut(NutArray, Var) REAL(r8) :: NutArray(grd%im,grd%jm,grd%km) REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: eva REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: evc - INTEGER(I4) :: MyNEofs, offset + INTEGER(I4) :: MyNEofs, MyNEofs_nit, offset CHARACTER :: Var NutArray(:,:,:) = 0.0 @@ -52,21 +52,26 @@ subroutine veof_nut(NutArray, Var) if((drv%nut .eq.1) .and. (drv%multiv .eq. 0)) then my_km = grd%km if(Var .eq. 'N') then - MyNEofs = ros%neof_n3n + MyNEofs_nit = ros%neof_n3n + MyNEofs = MyNEofs_nit offset = ros%neof_chl + ! If density-nutrient EOFs are provided and dnc flag set + if (drv%dnc .eq. 1) then + MyNEofs = MyNEofs_nit + ros%neof_dnc + endif else MyNEofs = ros%neof_o2o - offset = ros%neof_chl + ros%neof_n3n + offset = ros%neof_chl + ros%neof_n3n + ros%neof_dnc endif else if((drv%nut .eq.0) .and. (drv%multiv .eq. 1)) then - if(Var .eq. 'N') then + if((Var .eq. 'N') .and. (drv%dnc .eq. 0)) 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,*) "Error! Only nitrate multivariate assimilation without dens-nit covariances implemented" write(drv%dia,*) "" - write(*,*) "Error! Only nitrate multivariate assimilation implemented! Aborting" + write(*,*) "Error! Only nitrate multivariate assimilation without dens-nit covariances implemented! Aborting" write(*,*) "" endif call MPI_Barrier(Var3DCommunicator, ierr) @@ -90,17 +95,20 @@ subroutine veof_nut(NutArray, Var) 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 + eva(:,1:MyNEofs_nit) = ros%eva_n3n + evc(:,1:my_km,MyNEofs_nit) = ros%evc_n3n + if (drv%dnc .eq. 1) then + eva(:,MyNEofs_nit+1:MyNEofs) = ros%eva_dnc + evc(:,1:my_km,MyNEofs_nit+1:MyNEofs) = ros%evc_dnc(:,1:my_km,1:ros%neof_dnc) + endif 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 + if((Var .eq. 'N') .and. (drv%dnc.eq.0)) then eva = ros%eva_multi evc(:,1:my_km,:) = ros%evc_multi(:,ros%kmchl+1:ros%kmchl+grd%km,:) endif @@ -127,5 +135,5 @@ subroutine veof_nut(NutArray, Var) enddo enddo -DEALLOCATE(eva,evc) + DEALLOCATE(eva,evc) end subroutine veof_nut diff --git a/ver_hor_nut.f90 b/ver_hor_nut.f90 index be99ae2..08fc484 100644 --- a/ver_hor_nut.f90 +++ b/ver_hor_nut.f90 @@ -57,7 +57,11 @@ subroutine ver_hor_nut(NutArray, NutArrayAd, Var) ! --- ! Vertical EOFs - call veof_nut(NutArray, Var) + if(Var .eq. 'D') then + call veof_dnc(NutArray, Var) + else + call veof_nut(NutArray, Var) + endif !return ! goto 103 !No Vh diff --git a/ver_hor_nut_ad.f90 b/ver_hor_nut_ad.f90 index 8e02161..4b508aa 100644 --- a/ver_hor_nut_ad.f90 +++ b/ver_hor_nut_ad.f90 @@ -28,7 +28,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) INTEGER(i4) :: iProc, ierr type(DoubleGrid), allocatable, dimension(:,:,:) :: SendBuf3D type(DoubleGrid), allocatable, dimension(:) :: RecBuf1D - REAL(r8), allocatable, dimension(:,:,:) :: DefBufChl, DefBufChlAd + REAL(r8), allocatable, dimension(:,:,:) :: DefBufNut, DefBufNutAd REAL(r8) :: NutArray(grd%im,grd%jm,grd%km), NutArrayAd(grd%im,grd%jm,grd%km) CHARACTER :: Var @@ -65,20 +65,20 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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(DefBufNut(GlobalRow, localCol, grd%km)) + ALLOCATE(DefBufNutAd(GlobalRow, localCol, grd%km)) do k=1,grd%km do j=1,grd%jm do i=1,grd%im - SendBuf3D(k,i,j)%chl = NutArray(i,j,k) + SendBuf3D(k,i,j)%nut = NutArray(i,j,k) end do end do end do do k=1,grd%km do j=1,grd%jm do i=1,grd%im - SendBuf3D(k,i,j)%chl_ad = NutArrayAd(i,j,k) + SendBuf3D(k,i,j)%nut_ad = NutArrayAd(i,j,k) end do end do end do @@ -93,7 +93,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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 - DefBufChl(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl + DefBufNut(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%nut end do end do end do @@ -104,7 +104,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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 - DefBufChlAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl_ad + DefBufNutAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%nut_ad end do end do end do @@ -113,10 +113,10 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) ! --- ! Scale by the scaling factor do k=1,grd%km - DefBufChlAd(:,:,k) = DefBufChlAd(:,:,k) * grd%scx(:,:,k) + DefBufNutAd(:,:,k) = DefBufNutAd(:,:,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, grd%km, grd%imax, grd%aex, grd%bex, DefBufNutAd, grd%inx, grd%imx) else ! NumProcI .eq. 1 ! --- @@ -133,12 +133,12 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) ! 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, grd%km, grd%imax, grd%aex, grd%bex, DefBufNut, grd%inx, grd%imx) ! --- ! Scale by the scaling factor do k=1,grd%km - DefBufChl(:,:,k) = DefBufChl(:,:,k) * grd%scx(:,:,k) + DefBufNut(:,:,k) = DefBufNut(:,:,k) * grd%scx(:,:,k) enddo ! Reordering data to send back @@ -149,14 +149,14 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) do k=1,grd%km do j=1,localCol do i=1,GlobalRow - SendBuf3D(k,j,i)%chl = DefBufChl(i,j,k) + SendBuf3D(k,j,i)%nut = DefBufNut(i,j,k) end do end do end do do k=1,grd%km do j=1,localCol do i=1,GlobalRow - SendBuf3D(k,j,i)%chl_ad = DefBufChlAd(i,j,k) + SendBuf3D(k,j,i)%nut_ad = DefBufNutAd(i,j,k) end do end do end do @@ -171,7 +171,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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 - NutArray(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl + NutArray(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%nut end do end do end do @@ -182,13 +182,13 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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 - NutArrayAd(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl_ad + NutArrayAd(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%nut_ad end do end do end do end do - DEALLOCATE(SendBuf3D, RecBuf1D, DefBufChl, DefBufChlAd) + DEALLOCATE(SendBuf3D, RecBuf1D, DefBufNut, DefBufNutAd) else ! NumProcI .eq. 1 call rcfl_x( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, NutArray, grd%inx, grd%imx) @@ -223,7 +223,12 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) ! 103 continue ! --- ! Vertical EOFs - if(drv%multiv.eq.0) & - call veof_nut_ad(NutArrayAd, Var) + if(drv%multiv.eq.0) then + if(Var .eq. 'D') then + call veof_dnc_ad(NutArray, Var) + else + call veof_nut_ad(NutArrayAd, Var) + endif + endif end subroutine ver_hor_nut_ad From 513fa9ff0206a1780378607aa505b77a37654a70 Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Tue, 10 Feb 2026 22:18:49 +0100 Subject: [PATCH 08/13] Upadates to correct bugs --- Makefile | 7 +- def_nml.f90 | 2 +- grd_str.f90 | 2 + make3dvarwq.sh | 17 +++ obs_dnc_ad.f90 | 1 + resid.f90 | 1 + v_densnut.f90 | 91 ---------------- veof_dnc.f90 | 3 +- veof_nut.f90 | 6 +- ver_hor_densnut.f90 | 252 -------------------------------------------- ver_hor_nut_ad.f90 | 16 +-- 11 files changed, 41 insertions(+), 357 deletions(-) create mode 100644 make3dvarwq.sh delete mode 100644 v_densnut.f90 delete mode 100644 ver_hor_densnut.f90 diff --git a/Makefile b/Makefile index adaf875..4c826fd 100644 --- a/Makefile +++ b/Makefile @@ -72,6 +72,7 @@ OBJS = \ rdeofs_chl.o\ rdeofs_n3n.o\ rdeofs_o2o.o\ + rdeofs_dnc.o\ rdeofs_multi.o\ rdrcorr.o\ mean_rdr.o\ @@ -86,20 +87,22 @@ OBJS = \ cnv_ctv.o\ ver_hor_chl.o\ ver_hor_nut.o\ - ver_hor_densnut.o\ rcfl_x.o\ rcfl_y.o\ veof_chl.o\ veof_nut.o\ - v_densnut.o\ + veof_dnc.o\ obsop.o\ + obs_dnc.o\ obs_arg.o\ + obs_dnc_ad.o\ resid.o\ res_inc.o\ obsop_ad.o\ obs_arg_ad.o\ veof_chl_ad.o\ veof_nut_ad.o\ + veof_dnc_ad.o\ veof_multiv_ad.o\ ver_hor_chl_ad.o\ ver_hor_nut_ad.o\ diff --git a/def_nml.f90 b/def_nml.f90 index 4c92e1e..06b5bcf 100644 --- a/def_nml.f90 +++ b/def_nml.f90 @@ -42,7 +42,7 @@ subroutine def_nml implicit none LOGICAL :: read_eof - INTEGER(i4) :: neof_chl, neof_n3n, neof_o2o, nreg, rcf_ntr + INTEGER(i4) :: neof_chl, neof_n3n, neof_o2o, neof_dnc, nreg, rcf_ntr INTEGER(i4) :: neof_multi, kmchl, kmnit INTEGER(i4) :: verbose REAL(r8) :: rcf_L, ctl_tol, ctl_per, rcf_efc diff --git a/grd_str.f90 b/grd_str.f90 index c23b2b1..4415fd5 100644 --- a/grd_str.f90 +++ b/grd_str.f90 @@ -55,6 +55,8 @@ MODULE grd_str REAL(r8), POINTER :: n3n_ad(:,:,:) ! n3n adjoint variable REAL(r8), POINTER :: o2o(:,:,:) ! o2o REAL(r8), POINTER :: o2o_ad(:,:,:) ! o2o adjoint variable + REAL(r8), POINTER :: dnc(:,:,:) ! density + REAL(r8), POINTER :: dnc_ad(:,:,:) ! density adjoint variable REAL(r8), POINTER :: dep(:) ! Depth diff --git a/make3dvarwq.sh b/make3dvarwq.sh new file mode 100644 index 0000000..8e939ae --- /dev/null +++ b/make3dvarwq.sh @@ -0,0 +1,17 @@ +export MODULEFILE=$PWD/../ModelBuild/ogstm/compilers/machine_modules/g100.intel +source $MODULEFILE + + + +OGSTM_ARCH=x86_64 +OGSTM_OS=LINUX +OGSTM_COMPILER=intel + + +DEBUG_OCEANVAR=.dbg + + +INC_FILE=${OGSTM_ARCH}.${OGSTM_OS}.${OGSTM_COMPILER}${DEBUG_OCEANVAR}.inc +cp $INC_FILE compiler.inc +make clean +gmake diff --git a/obs_dnc_ad.f90 b/obs_dnc_ad.f90 index c55a01b..0627b7e 100644 --- a/obs_dnc_ad.f90 +++ b/obs_dnc_ad.f90 @@ -38,6 +38,7 @@ subroutine obs_dnc_ad use filenames use drv_str use bio_str + use dnc_str implicit none diff --git a/resid.f90 b/resid.f90 index 3253024..07acabc 100644 --- a/resid.f90 +++ b/resid.f90 @@ -32,6 +32,7 @@ subroutine resid use set_knd use obs_str use drv_str + use dnc_str implicit none diff --git a/v_densnut.f90 b/v_densnut.f90 deleted file mode 100644 index e33a4c5..0000000 --- a/v_densnut.f90 +++ /dev/null @@ -1,91 +0,0 @@ -subroutine v_densnut(NutArray) - - - !--------------------------------------------------------------------------- - ! ! - ! 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 . ! - ! ! - !--------------------------------------------------------------------------- - - !----------------------------------------------------------------------- - ! ! - ! Apply density nutrient increments and unterpolare aon vertical grid ! - ! ! - ! Version 1: A. Teruzzi 2025 ! - !----------------------------------------------------------------------- - - - use set_knd - use grd_str - ! use eof_str - ! use obs_str - use mpi_str - ! use filenames - use drv_str - use bio_str - use dnc_str - - implicit none - - INTEGER(i4) :: i, j, k, kk, my_km - REAL(r8), DIMENSION(grd%jm,grd%km) :: slicevar - REAL(r8), DIMENSION(grd%im+1,grd%jm,grd%km) :: N3nExt - REAL(r8) :: NutArray(grd%im,grd%jm,grd%km) - REAL(8) :: nit_dinc - - my_km = grd%km - - ! call EXTEND_2D( grd%n3n_ad, grd%km, N3nExtended_3d ) - - dnc%k = 0 - N3nExt(:,:,:) = 0.0 - - do kk = 1,dnc%no - - i=dnc%ib(kk) - j=dnc%jb(kk) - k=dnc%kb(kk) - - if(dnc%flc(kk).eq.1) then - - dnc%k = dnc%k + 1 - nit_dinc = dnc%inc(dnc%k) * dnc%corr(dnc%k) * dnc%err(dnc%k) / dnc%std(dnc%k) - - N3nExt(i ,j ,k ) = N3nExt(i ,j ,k ) + dnc%pq1(kk) * nit_dinc - N3nExt(i+1,j ,k ) = N3nExt(i+1,j ,k ) + dnc%pq2(kk) * nit_dinc - N3nExt(i ,j+1,k ) = N3nExt(i ,j+1,k ) + dnc%pq3(kk) * nit_dinc - N3nExt(i+1,j+1,k ) = N3nExt(i+1,j+1,k ) + dnc%pq4(kk) * nit_dinc - N3nExt(i ,j ,k+1) = N3nExt(i ,j ,k+1) + dnc%pq5(kk) * nit_dinc - N3nExt(i+1,j ,k+1) = N3nExt(i+1,j ,k+1) + dnc%pq6(kk) * nit_dinc - N3nExt(i ,j+1,k+1) = N3nExt(i ,j+1,k+1) + dnc%pq7(kk) * nit_dinc - N3nExt(i+1,j+1,k+1) = N3nExt(i+1,j+1,k+1) + dnc%pq8(kk) * nit_dinc - - - endif - - enddo - -! we apply contribution in grd%variable - - slicevar(:,1:grd%km) = NutArray(1,:,:) - call ADD_PREVCORE_CONTRIB(N3nExt, grd%km, NutArray, slicevar) - ! call ADD_PREVCORE_CONTRIB(N3nExtended_3d, grd%km, grd%N3n_ad, grd%n3n_ad(1,:,:)) - - - -end subroutine v_densnut diff --git a/veof_dnc.f90 b/veof_dnc.f90 index 26aa304..e957fdf 100644 --- a/veof_dnc.f90 +++ b/veof_dnc.f90 @@ -37,7 +37,8 @@ subroutine veof_dnc implicit none - INTEGER(i4) :: i, j, k, l,n, my_km, MyNEofs, ierr + INTEGER(i4) :: i, j, k, l,n, ierr + INTEGER(i4) :: my_km, MyNEofs, offset REAL(r8), DIMENSION ( grd%im, grd%jm) :: egm REAL(r8), ALLOCATABLE, DIMENSION(:,:) :: eva REAL(r8), ALLOCATABLE, DIMENSION(:,:,:) :: evc diff --git a/veof_nut.f90 b/veof_nut.f90 index abff91b..cbd0dfd 100644 --- a/veof_nut.f90 +++ b/veof_nut.f90 @@ -97,8 +97,10 @@ subroutine veof_nut(NutArray, Var) 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(:,1:MyNEofs_nit) = ros%eva_n3n - evc(:,1:my_km,MyNEofs_nit) = ros%evc_n3n + if(MyNEofs_nit .gt. 0) then + eva(:,1:MyNEofs_nit) = ros%eva_n3n + evc(:,1:my_km,1:MyNEofs_nit) = ros%evc_n3n + endif if (drv%dnc .eq. 1) then eva(:,MyNEofs_nit+1:MyNEofs) = ros%eva_dnc evc(:,1:my_km,MyNEofs_nit+1:MyNEofs) = ros%evc_dnc(:,1:my_km,1:ros%neof_dnc) diff --git a/ver_hor_densnut.f90 b/ver_hor_densnut.f90 deleted file mode 100644 index 3598586..0000000 --- a/ver_hor_densnut.f90 +++ /dev/null @@ -1,252 +0,0 @@ -subroutine ver_hor_densnut(NutArray) - - !--------------------------------------------------------------------------- - ! ! - ! 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 . ! - ! ! - !--------------------------------------------------------------------------- - - !----------------------------------------------------------------------- - ! ! - ! Apply horizontal filter ! - ! ! - ! Version 1: S.Dobricic 2006 ! - ! Version 2: S.Dobricic 2007 ! - ! Symmetric calculation in presence of coastal boundaries ! - ! eta_ad, tem_ad, and sal_ad are here temporary arrays ! - ! Version 3: A. Teruzzi 2013 ! - ! Attenuation of correction near the cost where d<200m ! - !----------------------------------------------------------------------- - - - use set_knd - use grd_str -! use eof_str - use cns_str - use drv_str - use dnc_str -! use obs_str - use mpi_str - - implicit none - - INTEGER(i4) :: i,j,k, ione, l - INTEGER :: jp, SurfaceIndex, TmpOffset, LinearIndex - INTEGER(i4) :: iProc, ierr - type(DoubleGrid), allocatable, dimension(:,:,:) :: SendBuf3D - type(DoubleGrid), allocatable, dimension(:) :: RecBuf1D(:) - REAL(r8), allocatable, dimension(:,:,:) :: DefBufChl, DefBufChlAd - REAL(r8) :: NutArray(grd%im,grd%jm,grd%km), NutArrayAd(grd%im,grd%jm,grd%km) -! CHARACTER :: Var - - ione = 1 - - ! --- - ! Apply density-nutrient increments and interp on vertical grid - call v_densnut(NutArray) - !return - ! goto 103 !No Vh - - ! --- - ! Load temporary arrays - do k=1,grd%km - NutArrayAd(:,:,k) = NutArray(:,:,k) - enddo - - !********** APPLY RECURSIVE FILTERS ********** ! - ! --- - ! Transpose calculation in the presense of coastal boundaries - - ! --- - ! y direction - ! --- - ! Scale by the scaling factor - do k=1,grd%km - NutArrayAd(:,:,k) = NutArrayAd(:,:,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, NutArrayAd, grd%jnx, grd%jmx) - - ! --- - ! 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)) - - do k=1,grd%km - do j=1,grd%jm - do i=1,grd%im - SendBuf3D(k,i,j)%chl = NutArray(i,j,k) - end do - end do - end do - do k=1,grd%km - do j=1,grd%jm - do i=1,grd%im - SendBuf3D(k,i,j)%chl_ad = NutArrayAd(i,j,k) - end do - end do - end do - - call MPI_Alltoallv(SendBuf3D, SendCountX3D, SendDisplX3D, MyPair, & - RecBuf1D, RecCountX3D, RecDisplX3D, MyPair, Var3DCommunicator, ierr) - - SurfaceIndex = localCol*grd%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 - DefBufChl(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl - end do - end do - - end do - 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 - DefBufChlAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl_ad - end do - end do - end do - end do - - ! --- - ! Scale by the scaling factor - do k=1,grd%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) - - else - ! --- - ! Scale by the scaling factor - do k=1,grd%km - NutArrayAd(:,:,k) = NutArrayAd(:,:,k) * grd%scx(:,:,k) - enddo - - call rcfl_x_ad( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, NutArrayAd, grd%inx, grd%imx) - - end if - - - - ! --- - ! 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) - - do k=1,grd%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)) - - do k=1,grd%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 j=1,localCol - do i=1,GlobalRow - SendBuf3D(k,j,i)%chl_ad = DefBufChlAd(i,j,k) - end do - end do - end do - - call MPI_Alltoallv(SendBuf3D, RecCountX3D, RecDisplX3D, MyPair, & - RecBuf1D, SendCountX3D, SendDisplX3D, MyPair, Var3DCommunicator, ierr) - - SurfaceIndex = grd%im*grd%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 - NutArray(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl - end do - end do - end do - 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 - NutArrayAd(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl_ad - end do - end do - end do - end do - - DEALLOCATE(SendBuf3D, RecBuf1D, DefBufChl, DefBufChlAd) - - else ! NumProcI .eq. 1 - - call rcfl_x( GlobalRow, localCol, grd%km, grd%imax, grd%aex, grd%bex, NutArray, grd%inx, grd%imx) - - do k=1,grd%km - NutArray(:,:,k) = NutArray(:,:,k) * grd%scx(:,:,k) - enddo - - end if - - ! --- - ! y direction - ! Apply recursive filter in y direction - call rcfl_y( localRow, GlobalCol, grd%km, grd%jmax, grd%aey, grd%bey, NutArray, grd%jnx, grd%jmx) - - ! --- - ! Scale by the scaling factor - do k=1,grd%km - NutArray(:,:,k) = NutArray(:,:,k) * grd%scy(:,:,k) - enddo - - ! --- - ! Average - do k=1,grd%km - NutArray(:,:,k) = (NutArray(:,:,k) + NutArrayAd(:,:,k) ) * 0.5 - enddo - - ! --- - ! Scale for boundaries - do k=1,grd%km - NutArray(:,:,k) = NutArray(:,:,k) * grd%msk(:,:,k) - enddo - - ! 103 continue - -end subroutine ver_hor_densnut diff --git a/ver_hor_nut_ad.f90 b/ver_hor_nut_ad.f90 index 4b508aa..c4b0d13 100644 --- a/ver_hor_nut_ad.f90 +++ b/ver_hor_nut_ad.f90 @@ -71,14 +71,14 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) do k=1,grd%km do j=1,grd%jm do i=1,grd%im - SendBuf3D(k,i,j)%nut = NutArray(i,j,k) + SendBuf3D(k,i,j)%chl = NutArray(i,j,k) end do end do end do do k=1,grd%km do j=1,grd%jm do i=1,grd%im - SendBuf3D(k,i,j)%nut_ad = NutArrayAd(i,j,k) + SendBuf3D(k,i,j)%chl_ad = NutArrayAd(i,j,k) end do end do end do @@ -93,7 +93,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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 - DefBufNut(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%nut + DefBufNut(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl end do end do end do @@ -104,7 +104,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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 - DefBufNutAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%nut_ad + DefBufNutAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl_ad end do end do end do @@ -149,14 +149,14 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) do k=1,grd%km do j=1,localCol do i=1,GlobalRow - SendBuf3D(k,j,i)%nut = DefBufNut(i,j,k) + SendBuf3D(k,j,i)%chl = DefBufNut(i,j,k) end do end do end do do k=1,grd%km do j=1,localCol do i=1,GlobalRow - SendBuf3D(k,j,i)%nut_ad = DefBufNutAd(i,j,k) + SendBuf3D(k,j,i)%chl_ad = DefBufNutAd(i,j,k) end do end do end do @@ -171,7 +171,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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 - NutArray(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%nut + NutArray(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl end do end do end do @@ -182,7 +182,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var) 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 - NutArrayAd(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%nut_ad + NutArrayAd(i, j + TmpOffset,k) = RecBuf1D(k + LinearIndex)%chl_ad end do end do end do From b014760d8429e17f28c3aa21780dbef017b7aa09 Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Thu, 12 Feb 2026 12:44:26 +0100 Subject: [PATCH 09/13] Small corrections on IFs and a filename --- filename_mod.f90 | 2 +- sav_itr.f90 | 5 ++++- wrt_dia.f90 | 37 +++++-------------------------------- wrt_nut_stat.f90 | 2 +- 4 files changed, 11 insertions(+), 35 deletions(-) diff --git a/filename_mod.f90 b/filename_mod.f90 index cf6bcb6..4a38e5f 100644 --- a/filename_mod.f90 +++ b/filename_mod.f90 @@ -43,7 +43,7 @@ SUBROUTINE SETFILENAMES RCORR_FILE = 'chl_rad_corr.nc' ARGO_FILE = 'arg_mis.dat' ANIS_FILE = 'gradsal.nc' - INCR_FILE = 'density_increments.nc' + INCR_FILE = 'density_increments.dat' END SUBROUTINE SETFILENAMES diff --git a/sav_itr.f90 b/sav_itr.f90 index 576c780..0b70097 100644 --- a/sav_itr.f90 +++ b/sav_itr.f90 @@ -114,6 +114,9 @@ subroutine sav_itr if(bio%O2o .eq. 1) then DEALLOCATE( ros%evc_o2o, ros%eva_o2o ) endif + if(drv%dnc .eq. 1) then + DEALLOCATE( ros%evc_dnc, ros%eva_dnc ) + endif endif if(drv%multiv.eq.1) then @@ -140,7 +143,7 @@ subroutine sav_itr endif if(drv%nut .eq. 1) then DEALLOCATE( bio%InitialNut) - if((bio%N3n.eq.1 .AND. bio%updateN1p.eq.1) .or. (drv%dnc .eq. 1)) DEALLOCATE(bio%covn3n_n1p) + if((bio%N3n.eq.1 .or. drv%dnc .eq. 1) .AND. (bio%updateN1p.eq.1)) DEALLOCATE(bio%covn3n_n1p) if(drv%chl_assim .eq. 0) then DEALLOCATE( bio%cquot, bio%pquot) DEALLOCATE( bio%InitialChl) !used in cp_chl_stat diff --git a/wrt_dia.f90 b/wrt_dia.f90 index 273da89..ad43c66 100644 --- a/wrt_dia.f90 +++ b/wrt_dia.f90 @@ -105,17 +105,19 @@ subroutine wrt_dia if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_att', status) endif - if(drv%dnc .eq. 1) then + if(drv%dnc .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) + if(bio%updateN1p .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 status = nf90mpi_enddef(ncid) if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', status) @@ -148,7 +150,7 @@ subroutine wrt_dia endif - if((drv%nut .eq. 1 .and. bio%n3n .eq. 1) .or. (drv%multiv .eq. 1)) then + if(( drv%nut .eq. 1 .and. ((bio%n3n .eq. 1) .or. (drv%dnc .eq. 1)) ) .or. (drv%multiv .eq. 1)) then do k=1,grd%km do j=1,grd%jm do i=1,grd%im @@ -165,7 +167,7 @@ subroutine wrt_dia endif if (bio%updateN1p .eq. 1) then - if((drv%nut .eq. 1 .and. bio%n3n .eq. 1).or.(drv%multiv.eq.1)) then + if(( drv%nut .eq. 1 .and. ((bio%n3n .eq. 1) .or. (drv%dnc .eq. 1)) ) .or. (drv%multiv .eq. 1)) then do k=1,grd%km do j=1,grd%jm do i=1,grd%im @@ -202,35 +204,6 @@ subroutine wrt_dia if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all o2o', status) endif - if(drv%dnc .eq. 1) then - do k=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) = REAL(grd%n3n(i,j,k), 4 ) - else - DumpMatrix(i,j,k) = 1.e20 - endif - enddo - enddo - enddo - status = nf90mpi_put_var_all(ncid,idn3n,DumpMatrix,MyStart,MyCount) - if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_var_all n3n', status) - - do k=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) = REAL(grd%n3n(i,j,k)*bio%covn3n_n1p(i,j,k), 4 ) - else - DumpMatrix(i,j,k) = 1.e20 - endif - enddo - enddo - enddo - 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 status = nf90mpi_close(ncid) if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_close', status) diff --git a/wrt_nut_stat.f90 b/wrt_nut_stat.f90 index cf82cc2..c2d668a 100644 --- a/wrt_nut_stat.f90 +++ b/wrt_nut_stat.f90 @@ -66,7 +66,7 @@ subroutine wrt_nut_stat call MPI_Abort(Var3DCommunicator,-1,ierr) endif - if((bio%updateN1p .eq. 1) .and. (NNVar.lt.2) .and. (drv%dnc .eq. 0)) then + if((bio%updateN1p .eq. 1) .and. (NNVar.lt.2) ) then write(*,*) "ERROR: Required phosphate update but NOT set in DA_params.f90" write(drv%dia,*) "ERROR: Required phosphate update but NOT set in DA_params.f90" call MPI_Barrier(Var3DCommunicator, ierr) From 437be94cedb984af50a5ace7445766b96d7b4285 Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Thu, 12 Feb 2026 12:50:06 +0100 Subject: [PATCH 10/13] Updated 3dvar namelist --- namelists/var_3d_nml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/namelists/var_3d_nml b/namelists/var_3d_nml index a2cbc7d..73ad2ec 100644 --- a/namelists/var_3d_nml +++ b/namelists/var_3d_nml @@ -28,6 +28,7 @@ ! 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_dnc - Number of vertical EOFs for density-nitrate ! neof_multi - Number of multivariate vertical EOFs ! kmchl - chl Nlevels in multi EOFs ! kmnit - nit Nlevels in multi EOFs @@ -43,6 +44,7 @@ neof_chl = 4 neof_n3n = 4 neof_o2o = 4 + neof_dnc = 10 neof_multi = 26 kmchl = 26 kmnit = 40 From 14f67d413346b0e2dcef832e81183441c6d1a617 Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Thu, 12 Feb 2026 14:07:30 +0100 Subject: [PATCH 11/13] onsistency checks with the help of chatGPT (see .md file) --- CODE_REVIEW_SESSION.md | 363 +++++++++++++++++++++++++++++++++++++++++ clean_mem.f90 | 1 + obs_arg.f90 | 2 +- obs_arg_ad.f90 | 2 +- readGrid.f90 | 12 +- sav_itr.f90 | 9 +- wrt_dia.f90 | 17 +- 7 files changed, 378 insertions(+), 28 deletions(-) create mode 100755 CODE_REVIEW_SESSION.md diff --git a/CODE_REVIEW_SESSION.md b/CODE_REVIEW_SESSION.md new file mode 100755 index 0000000..23175e0 --- /dev/null +++ b/CODE_REVIEW_SESSION.md @@ -0,0 +1,363 @@ +# Code Consistency Review Session - 3dVarBio with Density-Nutrient Coupling + +**Date**: February 12, 2026 + +--- + +## Overview + +This document contains the complete chat session reviewing the consistency of a modified version of the 3dVarBio marine data assimilation code that now uses density increments (from independent physical assimilation) to update nitrate (n3n). + +--- + +## User Context + +### Initial Description + +3dvarbio is a code to perform marine data assimilation for biogeochemistry. + +### Modification Description + +The user modified the code to consider increments on nitrate (usually named n3n in the code) based on density increments, which are resulting from an independent physical assimilation. + +### Previous Version Architecture + +The original version was designed to assimilate: +- **Surface chlorophyll** from satellite observations (sat) +- **Vertical profiles** of: + - Chlorophyll (chl) + - Oxygen (o2o/oxy) + - Nitrate (n3n/nit) + +In the new version, **nitrate increments are derived from density increments** from a separate physical assimilation, rather than being directly assimilated. + +--- + +## Code Consistency Review Request + +The user requested a consistency check focusing on three specific areas: + +### 1. Missing `drv%dnc` Conditionals + +**Question**: Do you see any missing if on `drv%dnc`, which is supposed to activate the density-nutrient assimilation? + +### 2. Observational Vector Consistency + +**Question**: `obs%amo` has to be consistent. This means that `obs%inc` and `obs%res` should be calculated on the correct variables (e.g., oxygen `obs%res` to be subtracted to the oxygen `obs%inc`). + +### 3. EOF Weight Ordering + +**Question**: EOFs are used to represent covariances on the vertical for nitrate, oxygen and chlorophyll but another set of eofs is now used by the code to consider multivariate vertical eofs to account for density-nutrient covariance. The elements of `ctl%x_c` should be calculated for the correct EOFs. They provide the weight for each eof. The user wants to be sure that `ctl%x_c` elements provide the weights in the correct order: +- First for chlorophyll EOFs +- Then for nitrate EOFs +- Then for multivariate density-nutrient EOFs +- Finally for oxygen EOFs + +--- + +## Code Review Findings + +### 1. Missing `drv%dnc` Conditionals - **CRITICAL ISSUES FOUND** + +#### [wrt_dia.f90](wrt_dia.f90) - **CRITICAL ISSUES** + +**Issue 1 - Lines 87-91**: Variable definition is missing `drv%dnc` check +```fortran +if((drv%nut .eq. 1 .and. bio%n3n .eq. 1) .or. (drv%multiv .eq. 1)) then +``` +Should be: +```fortran +if(((drv%nut .eq. 1 .and. bio%n3n .eq. 1) .or. (drv%dnc .eq. 1)) .or. (drv%multiv .eq. 1)) then +``` + +**Issue 2 - Lines 108-115**: Duplicate netCDF variable definition +The code defines 'n3n' twice when both conditions are true: +- Lines 87-91: If `(drv%nut .eq. 1 .and. bio%n3n .eq. 1)` +- Lines 108-115: If `drv%dnc .eq. 1` + +This will cause a netCDF error ("Variable already exists"). Lines 108-115 should be removed, as they become redundant with the corrected line 87. + +#### [obs_arg.f90](obs_arg.f90) - **LOGIC ISSUE** + +**Issue**: Line 55 doesn't check `drv%dnc` +```fortran +if ((drv%nut.eq.1 .and. bio%N3n.eq.1) .or. (drv%multiv.eq.1)) then +``` +Should be: +```fortran +if ((drv%nut.eq.1 .and. (bio%N3n.eq.1 .or. drv%dnc.eq.1)) .or. (drv%multiv.eq.1)) then +``` +Without this, `grd%n3n` will not be extended when `drv%dnc .eq. 1`, causing undefined behavior in the observational operator. + +#### [obs_arg_ad.f90](obs_arg_ad.f90) - **LOGIC ISSUE** + +**Issue**: Line 59 doesn't check `drv%dnc` for adjoint operation +```fortran +if ((drv%nut.eq.1 .and. bio%N3n.eq.1 ) .or. (drv%multiv.eq.1)) then +``` +Should be: +```fortran +if ((drv%nut.eq.1 .and. (bio%N3n.eq.1 .or. drv%dnc.eq.1)) .or. (drv%multiv.eq.1)) then +``` +Without this, the adjoint of the observational operator will not propagate gradients correctly. + +#### [readGrid.f90](readGrid.f90) - **CRITICAL DOUBLE ALLOCATION BUG** + +**Issue**: Lines 146-159 allocate `grd%n3n` and `grd%n3n_ad` twice +```fortran +if(bio%N3n .eq. 1) then + ALLOCATE ( grd%n3n(grd%im,grd%jm,grd%km) ) + ALLOCATE ( grd%n3n_ad(grd%im,grd%jm,grd%km) ) +endif +if(drv%dnc .eq. 1) then + ALLOCATE ( grd%n3n(grd%im,grd%jm,grd%km) ) ! DUPLICATE! + ALLOCATE ( grd%n3n_ad(grd%im,grd%jm,grd%km) ) ! DUPLICATE! + ... +endif +``` + +If both `bio%N3n .eq. 1` AND `drv%dnc .eq. 1`, the arrays are allocated twice, causing a fatal error. + +**Solution**: Modify line 154 to: +```fortran +if(drv%dnc .eq. 1 .and. bio%N3n .eq. 0) then +``` + +#### [sav_itr.f90](sav_itr.f90) - **CRITICAL DOUBLE DEALLOCATION BUG** + +**Issue**: Lines 78-79 and 95-96 deallocate `grd%n3n` twice +```fortran +if(bio%n3n .eq. 1) then + DEALLOCATE( grd%n3n) + DEALLOCATE( grd%n3n_ad) +endif +... +if(drv%dnc .eq. 1) then + DEALLOCATE( grd%n3n) ! DUPLICATE! +endif +``` + +If both `bio%n3n .eq. 1` AND `drv%dnc .eq. 1`, the arrays are deallocated twice, causing a segmentation fault. + +**Solution**: Modify the deallocation at line 95 to: +```fortran +if(drv%dnc .eq. 1 .and. bio%n3n .eq. 0) then + DEALLOCATE( grd%n3n) +endif +``` + +--- + +### 2. Observational Vector Consistency - **✓ VERIFIED CORRECT** + +Reviewed [resid.f90](resid.f90): + +The calculation of `obs%amo = (obs%inc - obs%res) / obs%err` is correctly applied to all observation types: + +- **Lines 49-50**: ARGO observations - correctly pairs `arg%inc` and `arg%res` ✓ +- **Lines 62-63**: Satellite observations - correctly pairs `sat%inc` and `sat%res` ✓ +- **Lines 75-76**: **Density observations** - correctly pairs `dnc%inc` and `dnc%res` ✓ + +All residual calculations in [obs_vec.f90](obs_vec.f90) correctly separate observations by type: +- ARGO residuals → `arg%res` → `obs%res` +- Satellite residuals → `sat%res` → `obs%res` +- Density residuals → `dnc%res` → `obs%res` + +**Conclusion**: The observational vector is **consistent and correct**. Each observation type maintains its residual independently until combined into the unified observation vector. + +--- + +### 3. EOF Weight Ordering - **✓ VERIFIED CORRECT** + +The EOF weight ordering in `ctl%x_c` is properly maintained throughout the code. + +#### Confirmed Ordering (indices in grd%ro and ctl%x_c): + +1. **Chlorophyll EOFs** (indices 1 → `ros%neof_chl`) + - Source: [veof_chl.f90](veof_chl.f90), line 87 + - Uses offset = 0 + +2. **Nitrate EOFs** (indices `ros%neof_chl+1` → `ros%neof_chl+ros%neof_n3n`) + - Source: [veof_nut.f90](veof_nut.f90), line 57 + - Uses offset = `ros%neof_chl` + +3. **Density-Nutrient EOFs** (indices `ros%neof_chl+ros%neof_n3n+1` → `ros%neof_chl+ros%neof_n3n+ros%neof_dnc`) + - Source: [veof_dnc.f90](veof_dnc.f90), line 50 + - Uses offset = `ros%neof_chl + ros%neof_n3n` + +4. **Oxygen EOFs** (indices `ros%neof_chl+ros%neof_n3n+ros%neof_dnc+1` → end) + - Source: [veof_nut.f90](veof_nut.f90), line 64 + - Uses offset = `ros%neof_chl + ros%neof_n3n + ros%neof_dnc` + +#### Verification in [def_cov.f90](def_cov.f90:492): +```fortran +ros%neof = ros%neof_multi + ros%neof_chl + ros%neof_nut +ros%neof_nut = ros%neof_n3n + ros%neof_o2o + ros%neof_dnc +``` + +The total EOF count correctly accounts for all types and the order is consistently maintained across matrix allocations and transformations. + +**Conclusion**: The EOF weight ordering is **correct and consistent** ✓ + +--- + +## Summary of Recommendations + +### Critical Issues (Must Fix): + +1. **[readGrid.f90](readGrid.f90), line 154** - Fix double allocation + - Change: `if(drv%dnc .eq. 1) then` → `if(drv%dnc .eq. 1 .and. bio%N3n .eq. 0) then` + +2. **[sav_itr.f90](sav_itr.f90), line 95** - Fix double deallocation + - Change: `if(drv%dnc .eq. 1) then` → `if(drv%dnc .eq. 1 .and. bio%n3n .eq. 0) then` + +3. **[wrt_dia.f90](wrt_dia.f90), lines 87 and 108-115** - Fix netCDF variable definition + - Update line 87 to include `drv%dnc` condition + - Delete lines 108-115 (duplicate n3n definition) + +4. **[obs_arg.f90](obs_arg.f90), line 55** - Add `drv%dnc` check + - Update condition to include `drv%dnc` + +5. **[obs_arg_ad.f90](obs_arg_ad.f90), line 59** - Add `drv%dnc` check + - Update condition to include `drv%dnc` + +### Verified Correct: + +- ✓ Observational vector consistency (`obs%amo`, `obs%inc`, `obs%res`) +- ✓ EOF weight ordering in `ctl%x_c` + +--- + +## Files Analyzed + +- bio_conv_ad.f90 / bio_conv.f90 +- bio_mod_ad.f90 / bio_mod.f90 +- bio_str.f90 +- biovar.f90 +- clean_mem.f90 +- cnv_ctv.f90 / cnv_ctv_ad.f90 +- cnv_inn.f90 +- costf.f90 +- def_cov.f90 +- def_nml.f90 / def_nml_multi.f90 +- **dnc_str.f90** +- **eof_str.f90** +- get_obs.f90 +- **get_densincr_arg.f90** +- get_obs_arg.f90 / get_obs_sat.f90 +- ini_cfn.f90 +- **obs_arg.f90 / obs_arg_ad.f90** ← Issues found +- **obs_dnc.f90 / obs_dnc_ad.f90** +- obs_vec.f90 +- **obsop.f90 / obsop_ad.f90** +- **readGrid.f90** ← Issues found +- **resid.f90** ← Verified correct +- **sav_itr.f90** ← Issues found +- veof_chl.f90 / veof_chl_ad.f90 +- **veof_dnc.f90 / veof_dnc_ad.f90** +- veof_nut.f90 / veof_nut_ad.f90 +- veof_multiv_ad.f90 +- **wrt_dia.f90** ← Issues found +- wrt_chl_stat.f90 +- wrt_nut_stat.f90 +- wrt_o2o_stat.f90 + +--- + +## Follow-up Review: Verification of Modifications + +After the user made modifications to address the critical issues, a second review was conducted. + +### Request + +The user asked for two additional checks: +1. **Check again `drv%dnc` occurrences** - Verify that all modifications were applied correctly +2. **Check allocation and deallocation consistency**: + - Consistency of allocated dimensions + - Ensure all needed deallocations are prescribed + +### Findings + +#### 1. `drv%dnc` Conditionals - **VERIFICATION COMPLETE** ✓ + +All critical `drv%dnc` checks have been properly added and verified in the following locations: + +- [wrt_dia.f90](wrt_dia.f90:87,94,140,157) - ✓ Now correctly includes `drv%dnc` +- [obs_arg.f90](obs_arg.f90:55) - ✓ Fixed to include `drv%dnc` check +- [obs_arg_ad.f90](obs_arg_ad.f90:58) - ✓ Fixed for adjoint operation +- [readGrid.f90](readGrid.f90:147,150) - ✓ Fixed (allocation properly conditional) +- [sav_itr.f90](sav_itr.f90:77,80,118,147) - ✓ Fixed (deallocation properly conditional) + +All conditional checks now properly distinguish between direct nitrate assimilation (`bio%N3n .eq. 1`) and density-derived nitrate updates (`drv%dnc .eq. 1`). + +#### 2. Allocation/Deallocation Consistency + +##### **Allocation Dimensions - ALL CORRECT** ✓ + +All allocations are properly dimensioned and consistent: +- `grd%n3n` and `grd%n3n_ad` (drv%multiv=0): dimension `(grd%im, grd%jm, grd%km)` ✓ +- `grd%n3n` and `grd%n3n_ad` (drv%multiv=1): dimension `(grd%im, grd%jm, ros%kmnit)` ✓ +- `grd%dnc` and `grd%dnc_ad`: dimension `(grd%im, grd%jm, grd%km)` ✓ +- `bio%phy` and `bio%phy_ad`: consistent dimensions across all paths ✓ +- `bio%cquot` and `bio%pquot`: allocated in [readChlStat.f90](readChlStat.f90:51-52) when needed ✓ +- `bio%InitialChl` and `bio%InitialNut`: allocated in readstat routines, deallocated in [sav_itr.f90](sav_itr.f90) ✓ + +##### **Deallocation Issues - ACTION REQUIRED** ⚠️ + +**Issue: Missing deallocation of `dnc%inc` in [clean_mem.f90](clean_mem.f90)** + +In [get_densincr_arg.f90](get_densincr_arg.f90:120), `dnc%inc` is allocated: +```fortran +ALLOCATE ( dnc%inc(dnc%no)) +``` + +But in [clean_mem.f90](clean_mem.f90:73-81), it is **NOT deallocated**. The current code shows: +```fortran +if(drv%dnc .eq. 1) then + DEALLOCATE ( dnc%flc) + DEALLOCATE ( dnc%res) + DEALLOCATE ( dnc%err) + DEALLOCATE ( dnc%ib, dnc%jb, dnc%kb) + DEALLOCATE ( dnc%pq1, dnc%pq2, dnc%pq3, dnc%pq4) + DEALLOCATE ( dnc%pq5, dnc%pq6, dnc%pq7, dnc%pq8) +endif +``` + +**Note**: `dnc%flg` is deallocated in [get_densincr_arg.f90](get_densincr_arg.f90:387), so it does NOT need to be deallocated in [clean_mem.f90](clean_mem.f90). +**BUT **: initially suggested to be done by chatGPT + +**Fix**: Add `DEALLOCATE ( dnc%inc)` in [clean_mem.f90](clean_mem.f90). + +##### **Deallocation Consistency - COMPLETE** ✓ (except above) + +All other allocations have corresponding deallocations: +- Grid data: ✓ +- Observations (argo, sat): ✓ +- Covariance structures (eofs): ✓ +- Bio structures: ✓ +- Control structures: ✓ +- Temporary work arrays: ✓ + +### Summary of Follow-up Review + +**Status**: ✓ Mostly Complete with Minor Fix Needed + +**Action Items**: +1. **[clean_mem.f90](clean_mem.f90:73-81)** - Add one missing deallocation: + ```fortran + if(drv%dnc .eq. 1) then + DEALLOCATE ( dnc%flc) + DEALLOCATE ( dnc%inc) ! ← ADD THIS + DEALLOCATE ( dnc%res) + DEALLOCATE ( dnc%err) + DEALLOCATE ( dnc%ib, dnc%jb, dnc%kb) + DEALLOCATE ( dnc%pq1, dnc%pq2, dnc%pq3, dnc%pq4) + DEALLOCATE ( dnc%pq5, dnc%pq6, dnc%pq7, dnc%pq8) + endif + ``` + +**Overall Assessment**: +- All `drv%dnc` conditionals are properly implemented and verified +- Allocation/deallocation is consistent and properly dimensioned +- Only one missing deallocation (`dnc%inc`) needs to be added to complete the implementation + diff --git a/clean_mem.f90 b/clean_mem.f90 index e1a0d56..0d19740 100644 --- a/clean_mem.f90 +++ b/clean_mem.f90 @@ -73,6 +73,7 @@ subroutine clean_mem if(drv%dnc .eq. 1) then DEALLOCATE ( dnc%flc) DEALLOCATE ( dnc%res) + DEALLOCATE ( dnc%inc) ! DEALLOCATE ( dnc%corr) DEALLOCATE ( dnc%err) ! DEALLOCATE ( dnc%std) diff --git a/obs_arg.f90 b/obs_arg.f90 index ace2e80..15dba81 100644 --- a/obs_arg.f90 +++ b/obs_arg.f90 @@ -52,7 +52,7 @@ subroutine obs_arg 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 + if ((drv%nut.eq.1 .and. ((bio%N3n.eq.1) .or. (drv%dnc.eq.1))) .or. (drv%multiv.eq.1)) then condn = 1 call EXTEND_2D( grd%n3n, grd%km, N3nExtended_3d ) endif diff --git a/obs_arg_ad.f90 b/obs_arg_ad.f90 index db4b78d..b449e53 100644 --- a/obs_arg_ad.f90 +++ b/obs_arg_ad.f90 @@ -55,7 +55,7 @@ subroutine obs_arg_ad 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 + if ((drv%nut.eq.1 .and. ((bio%N3n.eq.1) .or. (drv%dnc.eq.1)) ) .or. (drv%multiv.eq.1)) then call EXTEND_2D( grd%n3n_ad, grd%km, N3nExtended_3d ) condn = 1 endif diff --git a/readGrid.f90 b/readGrid.f90 index 2290f4c..d4246a6 100644 --- a/readGrid.f90 +++ b/readGrid.f90 @@ -144,20 +144,18 @@ subroutine readGrid 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(drv%nut .eq. 1) then - if(bio%N3n .eq. 1) then + if((bio%N3n .eq. 1) .or (drv%dnc .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%dnc .eq. 1) then + ALLOCATE ( grd%dnc(grd%im,grd%jm,grd%km) ) ; grd%dnc = huge(grd%dnc(1,1,1)) + ALLOCATE ( grd%dnc_ad(grd%im,grd%jm,grd%km) ) ; grd%dnc_ad = huge(grd%dnc_ad(1,1,1)) + endif 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 - if(drv%dnc .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)) - ALLOCATE ( grd%dnc(grd%im,grd%jm,grd%km) ) ; grd%dnc = huge(grd%dnc(1,1,1)) - ALLOCATE ( grd%dnc_ad(grd%im,grd%jm,grd%km) ) ; grd%dnc_ad = huge(grd%dnc_ad(1,1,1)) - endif endif else if(drv%multiv .eq. 1) then diff --git a/sav_itr.f90 b/sav_itr.f90 index 0b70097..7034701 100644 --- a/sav_itr.f90 +++ b/sav_itr.f90 @@ -74,9 +74,13 @@ subroutine sav_itr DEALLOCATE( grd%chl_ad) endif if(drv%nut .eq. 1) then - if(bio%n3n .eq. 1) then + if((bio%n3n .eq. 1) .or. (drv%dnc .eq. 1)) then DEALLOCATE( grd%n3n) DEALLOCATE( grd%n3n_ad) + if(drv%dnc .eq. 1) then + DEALLOCATE( grd%dnc) + DEALLOCATE( grd%dnc_ad) + endif endif if(bio%o2o .eq. 1) then DEALLOCATE( grd%o2o) @@ -93,9 +97,6 @@ subroutine sav_itr DEALLOCATE( grd%n3n_ad) endif - if(drv%dnc .eq. 1) then - DEALLOCATE( grd%n3n) - endif ! Observational vector DEALLOCATE( obs%inc, obs%amo, obs%res) diff --git a/wrt_dia.f90 b/wrt_dia.f90 index ad43c66..3131c2e 100644 --- a/wrt_dia.f90 +++ b/wrt_dia.f90 @@ -84,14 +84,14 @@ subroutine wrt_dia 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) .or. (drv%multiv .eq. 1)) then + if((drv%nut .eq. 1 .and. ((bio%n3n .eq. 1) .or (drv%dnc .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(bio%updateN1p .eq. 1) then - if((drv%nut .eq. 1 .and. bio%n3n .eq. 1) .or. (drv%multiv .eq. 1)) then + if((drv%nut .eq. 1 .and. ((bio%n3n .eq. 1) .or (drv%dnc .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) @@ -105,19 +105,6 @@ subroutine wrt_dia if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_put_att', status) endif - if(drv%dnc .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) - - if(bio%updateN1p .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 status = nf90mpi_enddef(ncid) if (status .ne. NF90_NOERR ) call handle_err('nf90mpi_def_var', status) From 2ab2769b447bbfd88061302282b15eee1ba2bae6 Mon Sep 17 00:00:00 2001 From: Anna Teruzzi Date: Fri, 13 Feb 2026 16:07:35 +0100 Subject: [PATCH 12/13] Debugging --- CODE_REVIEW_SESSION.md | 8 ++++++++ clean_mem.f90 | 2 -- def_nml.f90 | 2 +- readGrid.f90 | 2 +- wrt_dia.f90 | 4 ++-- 5 files changed, 12 insertions(+), 6 deletions(-) diff --git a/CODE_REVIEW_SESSION.md b/CODE_REVIEW_SESSION.md index 23175e0..88adfae 100755 --- a/CODE_REVIEW_SESSION.md +++ b/CODE_REVIEW_SESSION.md @@ -361,3 +361,11 @@ All other allocations have corresponding deallocations: - Allocation/deallocation is consistent and properly dimensioned - Only one missing deallocation (`dnc%inc`) needs to be added to complete the implementation + +# Check by anna running the code + +**Date**: February 13, 2026 + +The code crashed because 2 variables were deallocated twice: +dnc%res and dnc%err both in [obs_vec.f90] and in [clean_mem.f90] +(the correct one is in [obs_vec.f90]) diff --git a/clean_mem.f90 b/clean_mem.f90 index 0d19740..f173f65 100644 --- a/clean_mem.f90 +++ b/clean_mem.f90 @@ -72,10 +72,8 @@ subroutine clean_mem ! density increments if(drv%dnc .eq. 1) then DEALLOCATE ( dnc%flc) - DEALLOCATE ( dnc%res) DEALLOCATE ( dnc%inc) ! DEALLOCATE ( dnc%corr) - DEALLOCATE ( dnc%err) ! DEALLOCATE ( dnc%std) DEALLOCATE ( dnc%ib, dnc%jb, dnc%kb) DEALLOCATE ( dnc%pq1, dnc%pq2, dnc%pq3, dnc%pq4) diff --git a/def_nml.f90 b/def_nml.f90 index 06b5bcf..de404c4 100644 --- a/def_nml.f90 +++ b/def_nml.f90 @@ -48,7 +48,7 @@ subroutine def_nml REAL(r8) :: rcf_L, ctl_tol, ctl_per, rcf_efc 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 + NAMELIST /covlst/ neof_chl, neof_n3n, neof_o2o, neof_dnc, neof_multi, kmchl, kmnit, nreg, read_eof, rcf_ntr, rcf_L, rcf_efc ! ------------------------------------------------------------------- diff --git a/readGrid.f90 b/readGrid.f90 index d4246a6..76fee62 100644 --- a/readGrid.f90 +++ b/readGrid.f90 @@ -144,7 +144,7 @@ subroutine readGrid 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(drv%nut .eq. 1) then - if((bio%N3n .eq. 1) .or (drv%dnc .eq. 1)) then + if((bio%N3n .eq. 1) .or. (drv%dnc .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%dnc .eq. 1) then diff --git a/wrt_dia.f90 b/wrt_dia.f90 index 3131c2e..04b8215 100644 --- a/wrt_dia.f90 +++ b/wrt_dia.f90 @@ -84,14 +84,14 @@ subroutine wrt_dia 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) .or (drv%dnc .eq. 1))) .or. (drv%multiv .eq. 1)) then + if((drv%nut .eq. 1 .and. ((bio%n3n .eq. 1) .or. (drv%dnc .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(bio%updateN1p .eq. 1) then - if((drv%nut .eq. 1 .and. ((bio%n3n .eq. 1) .or (drv%dnc .eq. 1))) .or. (drv%multiv .eq. 1)) then + if((drv%nut .eq. 1 .and. ((bio%n3n .eq. 1) .or. (drv%dnc .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) From fe6e27f2e9678ad65a15a848c64f5b03eb159c22 Mon Sep 17 00:00:00 2001 From: AnnaT Date: Fri, 13 Feb 2026 16:10:52 +0100 Subject: [PATCH 13/13] Update CODE_REVIEW_SESSION.md with findings and recommendations Added a new section to the code review document detailing the findings and recommendations for the 3dVarBio code review session, including critical issues and verification of modifications. --- CODE_REVIEW_SESSION.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CODE_REVIEW_SESSION.md b/CODE_REVIEW_SESSION.md index 88adfae..149406f 100755 --- a/CODE_REVIEW_SESSION.md +++ b/CODE_REVIEW_SESSION.md @@ -367,5 +367,6 @@ All other allocations have corresponding deallocations: **Date**: February 13, 2026 The code crashed because 2 variables were deallocated twice: -dnc%res and dnc%err both in [obs_vec.f90] and in [clean_mem.f90] -(the correct one is in [obs_vec.f90]) +dnc%res and dnc%err both in [obs_vec.f90](obs_vec.f90) and in [clean_mem.f90](clean_mem.f90) +(the correct one is in [obs_vec.f90](obs_vec.f90)) +