diff --git a/CMakeLists.txt b/CMakeLists.txt
index c49e6d9..a024865 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -13,20 +13,13 @@ 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
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/CODE_REVIEW_SESSION.md b/CODE_REVIEW_SESSION.md
new file mode 100755
index 0000000..149406f
--- /dev/null
+++ b/CODE_REVIEW_SESSION.md
@@ -0,0 +1,372 @@
+# 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
+
+
+# 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](obs_vec.f90) and in [clean_mem.f90](clean_mem.f90)
+(the correct one is in [obs_vec.f90](obs_vec.f90))
+
diff --git a/Makefile b/Makefile
index a9b6bdc..4c826fd 100644
--- a/Makefile
+++ b/Makefile
@@ -60,45 +60,9 @@ OBJSTR = \
ctl_str.o\
rcfl_mod.o\
bio_str.o\
- mpi_str.o
+ mpi_str.o\
+ dnc_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\
@@ -108,6 +72,7 @@ OBJS = \
rdeofs_chl.o\
rdeofs_n3n.o\
rdeofs_o2o.o\
+ rdeofs_dnc.o\
rdeofs_multi.o\
rdrcorr.o\
mean_rdr.o\
@@ -115,6 +80,7 @@ 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\
@@ -125,14 +91,18 @@ OBJS = \
rcfl_y.o\
veof_chl.o\
veof_nut.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\
@@ -167,7 +137,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 94%
rename from oceanvar.f90
rename to biovar.f90
index 2185806..667b67f 100644
--- a/oceanvar.f90
+++ b/biovar.f90
@@ -1,4 +1,4 @@
-subroutine oceanvar
+subroutine biovar
!---------------------------------------------------------------------------
! !
@@ -112,10 +112,9 @@ subroutine oceanvar
! 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
- call cp_nut_stat
endif
endif
endif
@@ -125,18 +124,19 @@ subroutine oceanvar
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
call cp_nut_stat
- else if((bio%O2o .eq. 0) .and. (bio%N3n .eq. 1)) then
+ 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
+
call sav_itr
if(MyId .eq. 0) write(drv%dia,*) 'out of sav_itr '
@@ -146,4 +146,4 @@ subroutine oceanvar
!-----------------------------------------------------------------
if(MyId .eq. 0) close(drv%dia)
-end subroutine oceanvar
+end subroutine biovar
diff --git a/clean_mem.f90 b/clean_mem.f90
index a05255b..f173f65 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,17 @@ 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%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)
@@ -78,7 +90,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/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..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')
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 bf81477..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
@@ -531,8 +537,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 .or. drv%dnc .eq. 1) .AND. bio%updateN1p.eq.1) then
+ if(MyId .eq. 0) write(drv%dia,*) 'read nutcov'
call readNutCov
endif
if(drv%chl_assim.eq.0) then
@@ -545,6 +553,6 @@ subroutine def_cov
call readNutStat
if(bio%updateN1p.eq.1) &
call readNutCov
- endif
+ endif
end subroutine def_cov
diff --git a/def_nml.f90 b/def_nml.f90
index f99c065..de404c4 100644
--- a/def_nml.f90
+++ b/def_nml.f90
@@ -37,17 +37,18 @@ subroutine def_nml
use ctl_str
use mpi_str
use bio_str
+ use dnc_str
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
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
! -------------------------------------------------------------------
@@ -99,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
@@ -115,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 d5a0293..034be84 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, 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, uniformL, anisL
+ NAMELIST /params/ sat_obs, argo, dnc, uniformL, anisL
! -------------------------------------------------------------------
@@ -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,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,*) ' 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,*) '------------------------------------------------------------'
@@ -132,6 +133,7 @@ subroutine def_nml_multi
drv%sat_obs = sat_obs
drv%argo_obs = argo
+ drv%dnc = dnc
drv%uniformL = uniformL
drv%anisL = anisL
diff --git a/dnc_str.f90 b/dnc_str.f90
new file mode 100644
index 0000000..ead588c
--- /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 3e73d84..f44a4a3 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) :: 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)
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 dcf6e14..4a38e5f 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'
@@ -14,10 +15,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
@@ -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'
@@ -38,10 +40,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_densincr_arg.f90 b/get_densincr_arg.f90
new file mode 100644
index 0000000..6d02b2a
--- /dev/null
+++ b/get_densincr_arg.f90
@@ -0,0 +1,391 @@
+subroutine get_densincr_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 Density increments !
+ ! !
+ ! Version 1: S.Dobricic 2006 !
+ !-----------------------------------------------------------------------
+
+ use set_knd
+ use drv_str
+ use grd_str
+ use dnc_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, TmpLon, TmpLat
+ REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpDpt, TmpErr!, TmpStd
+ REAL(r8), ALLOCATABLE, DIMENSION(:) :: TmpInc!, TmpCorr
+ INTEGER(i4) :: GlobalDncNum, Counter, ierr
+ character(len=1024) :: filename
+
+ dnc%no = 0
+ dnc%nc = 0
+
+
+ ! ---
+ ! Allocate memory for observations
+ if(MyId .eq. 0) then
+ open(511,file=trim(INCR_FILE))
+ read(511,'(I4)') GlobalDncNum
+ write(drv%dia,*)'Number of density increments: ', GlobalDncNum
+ endif
+
+ call MPI_Bcast(GlobalDncNum, 1, MPI_INT, 0, Var3DCommunicator, ierr)
+
+ if(GlobalDncNum .eq. 0)then
+ if(MyId .eq. 0) &
+ close(511)
+ return
+ endif
+
+ 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 density increments
+ do k=1,GlobalDncNum
+ read (511,*) &
+ TmpFlc(k), & !TmpPar(k), &
+ TmpLon(k), TmpLat(k), &
+ TmpDpt(k), &!TmpTim(k), &
+ TmpInc(k), &!TmpCorr(k), &
+ TmpErr(k) !, TmpStd(k)
+ end do
+ 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(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(drv%Verbose .eq. 1) &
+ print*, "MyId", MyId, "has",Counter,"Density increments"
+
+ dnc%no = Counter
+
+ 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%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))
+ ALLOCATE ( dnc%pq5(dnc%no), dnc%pq6(dnc%no), dnc%pq7(dnc%no), dnc%pq8(dnc%no))
+
+ 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
+ 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%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%ino(Counter) = TmpIno(k)
+ ! endif
+ endif
+
+ enddo
+
+
+ ! ---
+ ! Initialise quality flag
+ dnc%flg(:) = 1
+ dnc%rb(:) = 0
+
+ ! ---
+! Vertical interpolation parameters
+ 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( 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
+ enddo
+
+
+ ! ---
+ ! Count good observations
+ dnc%nc = 0
+ do k=1,dnc%no
+ 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.
+ 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
+ dnc%flc(:) = dnc%flg(:)
+
+ DEALLOCATE( TmpFlc)!, TmpPar)
+ DEALLOCATE( TmpLon, TmpLat)
+ DEALLOCATE( TmpDpt)!, TmpTim)
+ DEALLOCATE( TmpInc)!, TmpTim)
+! DEALLOCATE( TmpCorr)!, TmpTim)
+ DEALLOCATE( TmpErr)
+! DEALLOCATE( TMPStd)
+! DEALLOCATE( TmpIno)
+
+end subroutine get_densincr_arg
+
+
+
+subroutine int_par_dnc
+
+ !-----------------------------------------------------------------------
+ ! !
+ ! 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 dnc_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(dnc%no.gt.0) then
+
+ dnc%flc(:) = dnc%flg(:)
+
+ ! ---
+ ! Horizontal interpolation parameters
+ do k = 1,dnc%no
+ do j=1,grd%jm-1
+ do i=1,grd%im-1
+ 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 + (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 + (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 = (dnc%lat(k) - grd%lat(1,1)) / grd%dlt + 1.0
+ ! j1 = int(q1)
+ ! 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
+ dnc%ib(k) = i1
+ dnc%jb(k) = j1
+ dnc%pb(k) = (p1-i1)
+ dnc%qb(k) = (q1-j1)
+ else
+ dnc%flc(k) = 0
+ endif
+ enddo
+
+ ! ---
+ ! Undefine masked for multigrid
+ 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.) dnc%flc(k) = 0
+ endif
+ enddo
+
+ ! ---
+ ! Horizontal interpolation parameters for each masked grid
+ 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=dnc%ib(k)
+ p1=dnc%pb(k)
+ j1=dnc%jb(k)
+ q1=dnc%qb(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)
+ 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 )
+ 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)
+ 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 )
+ 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=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)
+ 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 )
+ 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)
+ 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 )
+ 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=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(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
+
+
+
+ ! ---
+ ! Count good observations
+ dnc%nc = 0
+ do k=1,dnc%no
+ if(dnc%flc(k).eq.1)then
+ dnc%nc = dnc%nc + 1
+ endif
+ enddo
+
+ endif
+
+ 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 density increments: ',dnc%nc_global
+ print*,'Good density increments: ',dnc%nc_global
+ end if
+
+ DEALLOCATE ( dnc%flg)
+ DEALLOCATE ( dnc%lon, dnc%lat, dnc%dpt) !, dnc%tim)
+ DEALLOCATE ( dnc%pb, dnc%qb, dnc%rb)
+
+end subroutine int_par_dnc
diff --git a/get_obs.f90 b/get_obs.f90
index ed087a3..84b6687 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%dnc .eq. 1) &
+ call get_densincr_arg
+
end subroutine get_obs
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/int_par.f90 b/int_par.f90
index c6a533c..a5007fd 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%dnc .eq. 1) &
+ call int_par_dnc
+
+ end subroutine int_par
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
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/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/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
/
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
diff --git a/obs_arg.f90 b/obs_arg.f90
index 39d95f8..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/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..0627b7e
--- /dev/null
+++ b/obs_dnc_ad.f90
@@ -0,0 +1,97 @@
+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
+ use dnc_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_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
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/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..76fee62 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))
@@ -143,9 +144,13 @@ 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))
@@ -167,7 +172,7 @@ subroutine readGrid
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
@@ -184,7 +189,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..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
@@ -64,4 +65,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/sav_itr.f90 b/sav_itr.f90
index 6df2aff..7034701 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
@@ -73,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)
@@ -91,6 +96,7 @@ subroutine sav_itr
DEALLOCATE( grd%n3n)
DEALLOCATE( grd%n3n_ad)
endif
+
! Observational vector
DEALLOCATE( obs%inc, obs%amo, obs%res)
@@ -109,8 +115,15 @@ 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
+ DEALLOCATE( ros%evc_multi, ros%eva_multi)
+ endif
+
! Control structure
DEALLOCATE( ctl%x_c, ctl%g_c)
@@ -131,7 +144,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 .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
@@ -147,6 +160,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/veof_dnc.f90 b/veof_dnc.f90
new file mode 100644
index 0000000..e957fdf
--- /dev/null
+++ b/veof_dnc.f90
@@ -0,0 +1,82 @@
+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, 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
+
+ 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..cbd0dfd 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,22 @@ 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
+ 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)
+ 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 +137,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..c4b0d13 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,8 +65,8 @@ 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
@@ -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)%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
- DefBufChlAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl_ad
+ DefBufNutAd(i + TmpOffset,j,k) = RecBuf1D(k + LinearIndex)%chl_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)%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)%chl_ad = DefBufChlAd(i,j,k)
+ SendBuf3D(k,j,i)%chl_ad = DefBufNutAd(i,j,k)
end do
end do
end do
@@ -188,7 +188,7 @@ subroutine ver_hor_nut_ad(NutArray, NutArrayAd, Var)
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
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..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%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,6 +105,7 @@ subroutine wrt_dia
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)
@@ -136,7 +137,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
@@ -153,7 +154,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
@@ -169,26 +170,27 @@ 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
+
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..c2d668a 100644
--- a/wrt_nut_stat.f90
+++ b/wrt_nut_stat.f90
@@ -59,14 +59,14 @@ 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)
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) ) 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