MODULE asmbgc !!=========================================================================== !! *** MODULE asmbgc *** !! Routines and declarations for biogeochemical assimilation !!=========================================================================== !! History : 3.6 ! 2018-02 (D. Ford) Adapted from existing modules !!--------------------------------------------------------------------------- !! 'key_asminc' : assimilation increment interface !! 'key_top' : passive tracers !! 'key_hadocc' : HadOCC model !! 'key_medusa' : MEDUSA model !! 'key_roam' : MEDUSA extras for carbonate cycle !! 'key_fabm' : ERSEM model coupled via FABM !! 'key_karaml' : Kara mixed layer depth !!--------------------------------------------------------------------------- !! asm_bgc_check_options : check bgc assimilation options !! asm_bgc_init_incs : initialise bgc increments !! asm_bgc_init_bkg : initialise bgc background !! asm_bgc_bal_wri : write out bgc balancing increments !! asm_bgc_bkg_wri : write out bgc background !! phyto2d_asm_inc : apply the ocean colour increments !! phyto3d_asm_inc : apply the 3D phytoplankton increments !! pco2_asm_inc : apply the pCO2/fCO2 increments !! ph_asm_inc : apply the pH increments !! bgc3d_asm_inc : apply the generic 3D BGC increments !!--------------------------------------------------------------------------- USE par_kind, ONLY: & ! kind parameters & wp USE par_oce, ONLY: & ! domain array sizes & jpi, & & jpj, & & jpk USE iom ! i/o USE oce, ONLY: & ! active tracer variables & tsn USE zdfmxl, ONLY : & ! mixed layer depth #if defined key_karaml & hmld_kara, & & ln_kara, & #endif & hmld, & & hmld_tref, & & hmlp, & & hmlpt USE asmpar, ONLY: & ! assimilation parameters & c_asmbkg, & & c_asmbal, & & nitbkg_r, & & nitavgbkg_r, & & nitdin_r, & & nitiaustr_r, & & nitiaufin_r USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA & asm_phyto2d_bal_medusa USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC & asm_phyto2d_bal_hadocc USE asmphyto2dbal_ersem, ONLY: & ! phyto2d balancing for ERSEM & asm_phyto2d_bal_ersem USE asmpco2bal, ONLY: & ! pCO2 balancing & asm_pco2_bal #if defined key_top USE trc, ONLY: & ! passive tracer variables & trn, & & trb, & & nittrc000 USE par_trc, ONLY: & ! passive tracer parameters & jptra #endif #if defined key_medusa USE sms_medusa ! MEDUSA parameters USE par_medusa ! MEDUSA parameters USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system & mocsy_interface USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion & f2pCO2 USE sms_medusa, ONLY: & ! MEDUSA diagnostic variables & pgrow_avg, & & ploss_avg, & & phyt_avg, & & mld_max #elif defined key_hadocc USE par_hadocc ! HadOCC parameters USE had_bgc_const ! HadOCC parameters USE trc, ONLY: & ! HadOCC diagnostic variables & pgrow_avg, & & ploss_avg, & & phyt_avg, & & mld_max, & & HADOCC_CHL USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio & cchl_p #elif defined key_fabm USE par_fabm ! FABM-ERSEM parameters USE trc, ONLY: & ! FABM-ERSEM diagnostic variables & pgrow_avg, & & ploss_avg, & & phyt_avg, & & mld_max #endif IMPLICIT NONE PRIVATE PUBLIC asm_bgc_check_options ! called by asm_inc_init in asminc.F90 PUBLIC asm_bgc_init_incs ! called by asm_inc_init in asminc.F90 PRIVATE asm_bgc_read_incs_2d ! called by asm_bgc_init_incs PRIVATE asm_bgc_read_incs_3d ! called by asm_bgc_init_incs PUBLIC asm_bgc_init_bkg ! called by asm_inc_init in asminc.F90 PUBLIC asm_bgc_bal_wri ! called by nemo_gcm in nemogcm.F90 PUBLIC asm_bgc_bkg_alloc ! called by asm_bkg_wri in asmbkg.F90 PUBLIC asm_bgc_bkg_tavg ! called by asm_bkg_wri in asmbkg.F90 PUBLIC asm_bgc_bkg_wri ! called by asm_bkg_wri in asmbkg.F90 PUBLIC phyto2d_asm_inc ! called by bgc_asm_inc in asminc.F90 PUBLIC phyto3d_asm_inc ! called by bgc_asm_inc in asminc.F90 PUBLIC pco2_asm_inc ! called by bgc_asm_inc in asminc.F90 PUBLIC ph_asm_inc ! called by bgc_asm_inc in asminc.F90 PUBLIC bgc3d_asm_inc ! called by bgc_asm_inc in asminc.F90 LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of balancing incs LOGICAL, PUBLIC :: ln_phytobal = .FALSE. !: No phytoplankton balancing LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No surface total log10(chlorophyll) increment LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom log10(chlorophyll) increment LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment LOGICAL, PUBLIC :: ln_slchlnaninc = .FALSE. !: No surface nanophyto log10(chlorophyll) increment LOGICAL, PUBLIC :: ln_slchlpicinc = .FALSE. !: No surface picophyto log10(chlorophyll) increment LOGICAL, PUBLIC :: ln_slchldininc = .FALSE. !: No surface dinoflag log10(chlorophyll) increment LOGICAL, PUBLIC :: ln_schltotinc = .FALSE. !: No surface total chlorophyll increment LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total log10(phyto C) increment LOGICAL, PUBLIC :: ln_slphydiainc = .FALSE. !: No surface diatom log10(phyto C) increment LOGICAL, PUBLIC :: ln_slphynoninc = .FALSE. !: No surface non-diatom log10(phyto C) increment LOGICAL, PUBLIC :: ln_spco2inc = .FALSE. !: No surface pCO2 increment LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE. !: No surface fCO2 increment LOGICAL, PUBLIC :: ln_plchltotinc = .FALSE. !: No profile total log10(chlorophyll) increment LOGICAL, PUBLIC :: ln_pchltotinc = .FALSE. !: No profile total chlorophyll increment LOGICAL, PUBLIC :: ln_pno3inc = .FALSE. !: No profile nitrate increment LOGICAL, PUBLIC :: ln_psi4inc = .FALSE. !: No profile silicate increment LOGICAL, PUBLIC :: ln_pdicinc = .FALSE. !: No profile dissolved inorganic carbon increment LOGICAL, PUBLIC :: ln_palkinc = .FALSE. !: No profile alkalinity increment LOGICAL, PUBLIC :: ln_pphinc = .FALSE. !: No profile pH increment LOGICAL, PUBLIC :: ln_po2inc = .FALSE. !: No profile oxygen increment INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc assimilation ! 1) hmld - Turbocline/mixing depth ! [W points] ! 2) hmlp - Density criterion ! (0.01 kg/m^3 change from 10m) ! [W points] ! 3) hmld_kara - Kara MLD ! [Interpolated] ! 4) hmld_tref - Temperature criterion ! (0.2 K change from surface) ! [T points] ! 5) hmlpt - Density criterion ! (0.01 kg/m^3 change from 10m) ! [T points] REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc ! <= 0 no maximum applied (switch turned off) ! > 0 absolute chl inc capped at this value REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc ! slchltot inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchldia_bkginc ! slchldia inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlnan_bkginc ! slchlnan inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlpic_bkginc ! slchlpic inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchldin_bkginc ! slchldin inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: schltot_bkginc ! schltot inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphytot_bkginc ! slphytot inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphydia_bkginc ! slphydia inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphynon_bkginc ! slphynon inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sfco2_bkginc ! sfco2 inc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc ! spco2 inc REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: plchltot_bkginc ! plchltot inc REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pchltot_bkginc ! pchltot inc REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pno3_bkginc ! pno3 inc REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: psi4_bkginc ! psi4 inc REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pdic_bkginc ! pdic inc REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: palk_bkginc ! palk inc REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pph_bkginc ! pph inc REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: po2_bkginc ! po2 inc REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto2d_balinc ! Balancing incs from ocean colour REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto3d_balinc ! Balancing incs from p(l)chltot REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc ! Balancing incs from spco2/sfco2 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ph_balinc ! Balancing incs from pph REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg ! Background phyto growth REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg ! Background phyto loss REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg ! Background phyto conc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg ! Background max MLD REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg ! Background tracer state REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: chl_bkg ! Background chl REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cchl_p_bkg ! Background c:chl #if defined key_fabm REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: totalk_bkg ! Background total alkalinity #endif REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: pgrow_avg_tavg ! Time average of pgrow_avg REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: ploss_avg_tavg ! Time average of ploss_avg REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: phyt_avg_tavg ! Time average of phyt_avg REAL(wp), DIMENSION(:,:,:,:), SAVE, ALLOCATABLE :: trn_tavg ! Time average of trn #if defined key_hadocc REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: HADOCC_CHL_tavg ! Time average of HADOCC_CHL REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: cchl_p_tavg ! Time average of cchl_p #elif defined key_fabm REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: totalk_tavg ! Time average of O3_TA #endif # include "domzgr_substitute.h90" CONTAINS SUBROUTINE asm_bgc_check_options !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_check_options *** !! !! ** Purpose : check validity of biogeochemical assimilation options !! !! ** Method : check settings of logicals !! !! ** Action : call ctl_stop/ctl_warn if required !! !! References : asm_inc_init !!------------------------------------------------------------------------ #if ! defined key_top || ( ! defined key_hadocc && ! defined key_medusa && ! defined key_fabm ) CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', & & ' but no compatible biogeochemical model is available' ) #endif #if defined key_hadocc IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slphydiainc .OR. & & ln_slphynoninc .OR. ln_psi4inc .OR. ln_pphinc .OR. ln_po2inc ) THEN CALL ctl_stop( ' Cannot assimilate PFTs, Si4, pH or O2 into HadOCC' ) ENDIF #endif #if defined key_medusa IF ( .NOT. ln_foam_medusa ) THEN CALL ctl_stop( ' MEDUSA data assimilation options not turned on: set ln_foam_medusa' ) ENDIF #endif #if defined key_fabm IF ( ln_pphinc ) THEN CALL ctl_stop( ' Cannot currently assimilate pH into FABM-ERSEM' ) ENDIF #endif IF ( ( ln_phytobal ).AND. & & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. & & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. & & ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. & & ( .NOT. ln_schltotinc ).AND.( .NOT. ln_slphytotinc ).AND. & & ( .NOT. ln_slphydiainc ).AND.( .NOT. ln_slphynoninc ) ) THEN CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', & & ' if not assimilating ocean colour,', & & ' so ln_phytobal will be set to .false.') ln_phytobal = .FALSE. ENDIF IF ( ( ln_balwri ).AND. & & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. & & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. & & ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. & & ( .NOT. ln_schltotinc ).AND. & & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. & & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. & & ( .NOT. ln_pchltotinc ).AND.( .NOT. ln_pphinc ).AND. & & ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN CALL ctl_warn( ' Balancing increments not being calculated', & & ' so ln_balwri will be set to .false.') ln_balwri = .FALSE. ENDIF IF ( ln_spco2inc .AND. ln_sfco2inc ) THEN CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' ) ENDIF IF ( ln_slchltotinc .AND. ln_schltotinc ) THEN CALL ctl_stop( ' Can only assimilate surface log10(chlorophyll) or chlorophyll, not both' ) ENDIF IF ( ln_plchltotinc .AND. ln_pchltotinc ) THEN CALL ctl_stop( ' Can only assimilate profile log10(chlorophyll) or chlorophyll, not both' ) ENDIF IF ( ( ln_slchltotinc .OR. ln_schltotinc ) .AND. & & ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slchlnaninc .OR. & & ln_slchlpicinc .OR. ln_slchldininc ) ) THEN CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' ) ENDIF IF ( ln_phytobal .AND. & & ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. & #if defined key_fabm & ( ( ln_slchldiainc .OR. ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc ) .AND. & & ( ( .NOT. ln_slchldiainc ) .OR. ( .NOT. ln_slchlnaninc ) .OR. & & ( .NOT. ln_slchlpicinc ) .OR. ( .NOT. ln_slchldininc ) ) ) .OR. & #endif & ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', & & ' unless assimilating all model PFTs,') ENDIF IF ( ln_slphytotinc .AND. ( ln_slphydiainc .OR. ln_slphynoninc ) ) THEN CALL ctl_stop( ' Can only assimilate total or PFT surface phytoplankton carbon, not both' ) ENDIF IF ( ln_phytobal .AND. & & ( ( ln_slphynoninc .AND.( .NOT. ln_slphydiainc ) ).OR. & & ( ln_slphydiainc .AND.( .NOT. ln_slphynoninc ) ) ) ) THEN CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', & & ' unless assimilating all model PFTs,') ENDIF END SUBROUTINE asm_bgc_check_options !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_init_incs( knum ) !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_init_incs *** !! !! ** Purpose : initialise bgc increments !! !! ** Method : allocate arrays and read increments from file !! !! ** Action : allocate arrays and read increments from file !! !! References : asm_inc_init !!------------------------------------------------------------------------ !! INTEGER, INTENT(in ) :: knum ! i/o unit of increments file !! !!------------------------------------------------------------------------ ! Allocate and read increments IF ( ln_slchltotinc ) THEN ALLOCATE( slchltot_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc ) ENDIF IF ( ln_slchldiainc ) THEN ALLOCATE( slchldia_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc ) ENDIF IF ( ln_slchlnoninc ) THEN ALLOCATE( slchlnon_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc ) ENDIF IF ( ln_slchlnaninc ) THEN ALLOCATE( slchlnan_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnan', slchlnan_bkginc ) ENDIF IF ( ln_slchlpicinc ) THEN ALLOCATE( slchlpic_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslchlpic', slchlpic_bkginc ) ENDIF IF ( ln_slchldininc ) THEN ALLOCATE( slchldin_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslchldin', slchldin_bkginc ) ENDIF IF ( ln_schltotinc ) THEN ALLOCATE( schltot_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc ) ENDIF IF ( ln_slphytotinc ) THEN ALLOCATE( slphytot_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc ) ENDIF IF ( ln_slphydiainc ) THEN ALLOCATE( slphydia_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc ) ENDIF IF ( ln_slphynoninc ) THEN ALLOCATE( slphynon_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinslphynon', slphynon_bkginc ) ENDIF IF ( ln_sfco2inc ) THEN ALLOCATE( sfco2_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinsfco2', sfco2_bkginc ) ENDIF IF ( ln_spco2inc ) THEN ALLOCATE( spco2_bkginc(jpi,jpj) ) CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc ) ENDIF IF ( ln_plchltotinc ) THEN ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) ) CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc ) ENDIF IF ( ln_pchltotinc ) THEN ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) ) CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc ) ENDIF IF ( ln_pno3inc ) THEN ALLOCATE( pno3_bkginc(jpi,jpj,jpk) ) CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc ) ENDIF IF ( ln_psi4inc ) THEN ALLOCATE( psi4_bkginc(jpi,jpj,jpk) ) CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc ) ENDIF IF ( ln_pdicinc ) THEN ALLOCATE( pdic_bkginc(jpi,jpj,jpk) ) CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc ) ENDIF IF ( ln_palkinc ) THEN ALLOCATE( palk_bkginc(jpi,jpj,jpk) ) CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc ) ENDIF IF ( ln_pphinc ) THEN ALLOCATE( pph_bkginc(jpi,jpj,jpk) ) CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc ) ENDIF IF ( ln_po2inc ) THEN ALLOCATE( po2_bkginc(jpi,jpj,jpk) ) CALL asm_bgc_read_incs_3d( knum, 'bckinpo2', po2_bkginc ) ENDIF ! Allocate balancing increments IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. & & ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & & ln_slphynoninc ) THEN #if defined key_top ALLOCATE( phyto2d_balinc(jpi,jpj,jpk,jptra) ) phyto2d_balinc(:,:,:,:) = 0.0 #else CALL ctl_stop( ' key_top must be set for balancing increments' ) #endif ENDIF IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN #if defined key_top ALLOCATE( phyto3d_balinc(jpi,jpj,jpk,jptra) ) phyto3d_balinc(:,:,:,:) = 0.0 #else CALL ctl_stop( ' key_top must be set for balancing increments' ) #endif ENDIF IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN #if defined key_top ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) ) pco2_balinc(:,:,:,:) = 0.0 #else CALL ctl_stop( ' key_top must be set for balancing increments' ) #endif ENDIF IF ( ln_pphinc ) THEN #if defined key_top ALLOCATE( ph_balinc(jpi,jpj,jpk,jptra) ) ph_balinc(:,:,:,:) = 0.0 #else CALL ctl_stop( ' key_top must be set for balancing increments' ) #endif ENDIF END SUBROUTINE asm_bgc_init_incs !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_read_incs_2d( knum, cd_bgcname, p_incs ) !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_init_incs *** !! !! ** Purpose : read 2d bgc increments !! !! ** Method : read increments from file !! !! ** Action : read increments from file !! !! References : asm_inc_init !!------------------------------------------------------------------------ !! INTEGER, INTENT(in ) :: knum ! i/o unit CHARACTER(LEN=*), INTENT(in ) :: cd_bgcname ! variable REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_incs ! increments !! !!------------------------------------------------------------------------ ! Initialise p_incs(:,:) = 0.0 ! read from file CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 ) ! Apply the masks p_incs(:,:) = p_incs(:,:) * tmask(:,:,1) ! Set missing increments to 0.0 rather than 1e+20 ! to allow for differences in masks WHERE( ABS( p_incs(:,:) ) > 1.0e+10 ) p_incs(:,:) = 0.0 END SUBROUTINE asm_bgc_read_incs_2d !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_read_incs_3d( knum, cd_bgcname, p_incs ) !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_init_incs *** !! !! ** Purpose : read 3d bgc increments !! !! ** Method : read increments from file !! !! ** Action : read increments from file !! !! References : asm_inc_init !!------------------------------------------------------------------------ !! INTEGER, INTENT(in ) :: knum ! i/o unit CHARACTER(LEN=*), INTENT(in ) :: cd_bgcname ! variable REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: p_incs ! increments !! !!------------------------------------------------------------------------ ! Initialise p_incs(:,:,:) = 0.0 ! read from file CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 ) ! Apply the masks p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:) ! Set missing increments to 0.0 rather than 1e+20 ! to allow for differences in masks WHERE( ABS( p_incs(:,:,:) ) > 1.0e+10 ) p_incs(:,:,:) = 0.0 END SUBROUTINE asm_bgc_read_incs_3d !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_init_bkg !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_init_bkg *** !! !! ** Purpose : initialise bgc background !! !! ** Method : allocate arrays and read background from file !! !! ** Action : allocate arrays and read background from file !! !! References : asm_inc_init !!------------------------------------------------------------------------ !! INTEGER :: inum ! i/o unit of background file INTEGER :: jt ! loop counter !! !!------------------------------------------------------------------------ #if defined key_top && ( defined key_hadocc || defined key_medusa || defined key_fabm ) IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. & & ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & & ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN ALLOCATE( pgrow_avg_bkg(jpi,jpj) ) ALLOCATE( ploss_avg_bkg(jpi,jpj) ) ALLOCATE( phyt_avg_bkg(jpi,jpj) ) ALLOCATE( mld_max_bkg(jpi,jpj) ) ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) pgrow_avg_bkg(:,:) = 0.0 ploss_avg_bkg(:,:) = 0.0 phyt_avg_bkg(:,:) = 0.0 mld_max_bkg(:,:) = 0.0 tracer_bkg(:,:,:,:) = 0.0 #if defined key_hadocc ALLOCATE( chl_bkg(jpi,jpj,jpk) ) ALLOCATE( cchl_p_bkg(jpi,jpj,jpk) ) chl_bkg(:,:,:) = 0.0 cchl_p_bkg(:,:,:) = 0.0 #endif !-------------------------------------------------------------------- ! Read background variables for phytoplankton assimilation ! Some only required if performing balancing !-------------------------------------------------------------------- CALL iom_open( c_asmbkg, inum ) #if defined key_hadocc CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl', chl_bkg ) CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg ) chl_bkg(:,:,:) = chl_bkg(:,:,:) * tmask(:,:,:) cchl_p_bkg(:,:,:) = cchl_p_bkg(:,:,:) * tmask(:,:,:) #elif defined key_medusa CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) ) #elif defined key_fabm CALL iom_get( inum, jpdom_autoglo, 'ersem_chl1', tracer_bkg(:,:,:,jp_fabm_chl1) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_chl2', tracer_bkg(:,:,:,jp_fabm_chl2) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_chl3', tracer_bkg(:,:,:,jp_fabm_chl3) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_chl4', tracer_bkg(:,:,:,jp_fabm_chl4) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p1c', tracer_bkg(:,:,:,jp_fabm_p1c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p1n', tracer_bkg(:,:,:,jp_fabm_p1n) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p1p', tracer_bkg(:,:,:,jp_fabm_p1p) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p1s', tracer_bkg(:,:,:,jp_fabm_p1s) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p2c', tracer_bkg(:,:,:,jp_fabm_p2c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p2n', tracer_bkg(:,:,:,jp_fabm_p2n) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p2p', tracer_bkg(:,:,:,jp_fabm_p2p) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p3c', tracer_bkg(:,:,:,jp_fabm_p3c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p3n', tracer_bkg(:,:,:,jp_fabm_p3n) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p3p', tracer_bkg(:,:,:,jp_fabm_p3p) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p4c', tracer_bkg(:,:,:,jp_fabm_p4c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p4n', tracer_bkg(:,:,:,jp_fabm_p4n) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_p4p', tracer_bkg(:,:,:,jp_fabm_p4p) ) #endif IF ( ln_phytobal ) THEN CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg ) CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg ) CALL iom_get( inum, jpdom_autoglo, 'phyt_avg', phyt_avg_bkg ) CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg ) pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1) ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1) phyt_avg_bkg(:,:) = phyt_avg_bkg(:,:) * tmask(:,:,1) mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) #if defined key_hadocc CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) ) CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) ) CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) ) CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) ) CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) #elif defined key_medusa CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) ) #elif defined key_fabm CALL iom_get( inum, jpdom_autoglo, 'ersem_z4c', tracer_bkg(:,:,:,jp_fabm_z4c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_z5c', tracer_bkg(:,:,:,jp_fabm_z5c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_z5n', tracer_bkg(:,:,:,jp_fabm_z5n) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_z5p', tracer_bkg(:,:,:,jp_fabm_z5p) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_z6c', tracer_bkg(:,:,:,jp_fabm_z6c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_z6n', tracer_bkg(:,:,:,jp_fabm_z6n) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_z6p', tracer_bkg(:,:,:,jp_fabm_z6p) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_n1p', tracer_bkg(:,:,:,jp_fabm_n1p) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_n3n', tracer_bkg(:,:,:,jp_fabm_n3n) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_n4n', tracer_bkg(:,:,:,jp_fabm_n4n) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_n5s', tracer_bkg(:,:,:,jp_fabm_n5s) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_o2o', tracer_bkg(:,:,:,jp_fabm_o2o) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) ) totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:) #endif ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN #if defined key_hadocc CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) #elif defined key_medusa CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) #elif defined key_fabm CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) ) totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:) #endif CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg ) mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) ENDIF CALL iom_close( inum ) DO jt = 1, jptra tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) END DO ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) ALLOCATE( mld_max_bkg(jpi,jpj) ) tracer_bkg(:,:,:,:) = 0.0 mld_max_bkg(:,:) = 0.0 CALL iom_open( c_asmbkg, inum ) #if defined key_hadocc CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) #elif defined key_medusa CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) #elif defined key_fabm CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) ) CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) ) totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:) #endif CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg ) CALL iom_close( inum ) DO jt = 1, jptra tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) END DO mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) ENDIF #else CALL ctl_stop( ' asm_bgc_init_bkg: key_top and a compatible biogeochemical model required' ) #endif END SUBROUTINE asm_bgc_init_bkg !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_bal_wri( kt ) !!------------------------------------------------------------------------ !! !! *** ROUTINE asm_bgc_bal_wri *** !! !! ** Purpose : Write to file the balancing increments !! calculated for biogeochemistry !! !! ** Method : Write to file the balancing increments !! calculated for biogeochemistry !! !! ** Action : !! !! References : !! !! History : !! ! 2014-08 (D. Ford) Initial version, based on asm_bkg_wri !!------------------------------------------------------------------------ !! * Arguments INTEGER, INTENT( IN ) :: kt ! Current time-step !! * Local declarations CHARACTER(LEN=50) :: cl_asmbal ! Filename (with extension) LOGICAL :: llok ! Check if file exists INTEGER :: inum ! File unit number REAL(wp) :: zdate ! Date !!------------------------------------------------------------------------ ! Set things up zdate = REAL( ndastp ) WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal ) ! Check if file exists INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok ) IF( .NOT. llok ) THEN IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// & & TRIM( c_asmbal ) // ' at timestep = ', kt ! Define the output file CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib) ! Write the information CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate ) IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. & & ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & & ln_slphynoninc ) THEN #if defined key_medusa CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) ) CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) ) CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) ) CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) ) CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) ) IF ( ln_phytobal ) THEN CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) ) CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) ) CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) ) CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) ) CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) ) CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) ) CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) ) CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) ) CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) ) CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) ) ENDIF #elif defined key_hadocc CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) ) IF ( ln_phytobal ) THEN CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) ) CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) ) CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) ) CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) ) ENDIF #elif defined key_fabm CALL iom_rstput( kt, kt, inum, 'phy2d_chl1', phyto2d_balinc(:,:,:,jp_fabm_chl1) ) CALL iom_rstput( kt, kt, inum, 'phy2d_chl2', phyto2d_balinc(:,:,:,jp_fabm_chl2) ) CALL iom_rstput( kt, kt, inum, 'phy2d_chl3', phyto2d_balinc(:,:,:,jp_fabm_chl3) ) CALL iom_rstput( kt, kt, inum, 'phy2d_chl4', phyto2d_balinc(:,:,:,jp_fabm_chl4) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p1c', phyto2d_balinc(:,:,:,jp_fabm_p1c) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p1n', phyto2d_balinc(:,:,:,jp_fabm_p1n) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p1p', phyto2d_balinc(:,:,:,jp_fabm_p1p) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p1s', phyto2d_balinc(:,:,:,jp_fabm_p1s) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p2c', phyto2d_balinc(:,:,:,jp_fabm_p2c) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p2n', phyto2d_balinc(:,:,:,jp_fabm_p2n) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p2p', phyto2d_balinc(:,:,:,jp_fabm_p2p) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p3c', phyto2d_balinc(:,:,:,jp_fabm_p3c) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p3n', phyto2d_balinc(:,:,:,jp_fabm_p3n) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p3p', phyto2d_balinc(:,:,:,jp_fabm_p3p) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p4c', phyto2d_balinc(:,:,:,jp_fabm_p4c) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p4n', phyto2d_balinc(:,:,:,jp_fabm_p4n) ) CALL iom_rstput( kt, kt, inum, 'phy2d_p4p', phyto2d_balinc(:,:,:,jp_fabm_p4p) ) IF ( ln_phytobal ) THEN CALL iom_rstput( kt, kt, inum, 'phy2d_z4c', phyto2d_balinc(:,:,:,jp_fabm_z4c) ) CALL iom_rstput( kt, kt, inum, 'phy2d_z5c', phyto2d_balinc(:,:,:,jp_fabm_z5c) ) CALL iom_rstput( kt, kt, inum, 'phy2d_z5n', phyto2d_balinc(:,:,:,jp_fabm_z5n) ) CALL iom_rstput( kt, kt, inum, 'phy2d_z5p', phyto2d_balinc(:,:,:,jp_fabm_z5p) ) CALL iom_rstput( kt, kt, inum, 'phy2d_z6c', phyto2d_balinc(:,:,:,jp_fabm_z6c) ) CALL iom_rstput( kt, kt, inum, 'phy2d_z6n', phyto2d_balinc(:,:,:,jp_fabm_z6n) ) CALL iom_rstput( kt, kt, inum, 'phy2d_z6p', phyto2d_balinc(:,:,:,jp_fabm_z6p) ) CALL iom_rstput( kt, kt, inum, 'phy2d_n1p', phyto2d_balinc(:,:,:,jp_fabm_n1p) ) CALL iom_rstput( kt, kt, inum, 'phy2d_n3n', phyto2d_balinc(:,:,:,jp_fabm_n3n) ) CALL iom_rstput( kt, kt, inum, 'phy2d_n4n', phyto2d_balinc(:,:,:,jp_fabm_n4n) ) CALL iom_rstput( kt, kt, inum, 'phy2d_n5s', phyto2d_balinc(:,:,:,jp_fabm_n5s) ) CALL iom_rstput( kt, kt, inum, 'phy2d_o2o', phyto2d_balinc(:,:,:,jp_fabm_o2o) ) CALL iom_rstput( kt, kt, inum, 'phy2d_o3c', phyto2d_balinc(:,:,:,jp_fabm_o3c) ) CALL iom_rstput( kt, kt, inum, 'phy2d_o3ba', phyto2d_balinc(:,:,:,jp_fabm_o3ba) ) ENDIF #endif ENDIF IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN #if defined key_medusa CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) ) CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) ) CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) ) CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) ) CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) ) #elif defined key_hadocc CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) ) #elif defined key_fabm CALL iom_rstput( kt, kt, inum, 'phy3d_chl1', phyto3d_balinc(:,:,:,jp_fabm_chl1) ) CALL iom_rstput( kt, kt, inum, 'phy3d_chl2', phyto3d_balinc(:,:,:,jp_fabm_chl2) ) CALL iom_rstput( kt, kt, inum, 'phy3d_chl3', phyto3d_balinc(:,:,:,jp_fabm_chl3) ) CALL iom_rstput( kt, kt, inum, 'phy3d_chl4', phyto3d_balinc(:,:,:,jp_fabm_chl4) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p1c', phyto3d_balinc(:,:,:,jp_fabm_p1c) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p1n', phyto3d_balinc(:,:,:,jp_fabm_p1n) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p1p', phyto3d_balinc(:,:,:,jp_fabm_p1p) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p1s', phyto3d_balinc(:,:,:,jp_fabm_p1s) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p2c', phyto3d_balinc(:,:,:,jp_fabm_p2c) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p2n', phyto3d_balinc(:,:,:,jp_fabm_p2n) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p2p', phyto3d_balinc(:,:,:,jp_fabm_p2p) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p3c', phyto3d_balinc(:,:,:,jp_fabm_p3c) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p3n', phyto3d_balinc(:,:,:,jp_fabm_p3n) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p3p', phyto3d_balinc(:,:,:,jp_fabm_p3p) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p4c', phyto3d_balinc(:,:,:,jp_fabm_p4c) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p4n', phyto3d_balinc(:,:,:,jp_fabm_p4n) ) CALL iom_rstput( kt, kt, inum, 'phy3d_p4p', phyto3d_balinc(:,:,:,jp_fabm_p4p) ) #endif ENDIF IF ( ln_spco2inc ) THEN #if defined key_medusa CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) ) CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) ) #elif defined key_hadocc CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) ) #elif defined key_fabm CALL iom_rstput( kt, kt, inum, 'pco2_o3c', pco2_balinc(:,:,:,jp_fabm_o3c) ) CALL iom_rstput( kt, kt, inum, 'pco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) ) #endif ELSE IF ( ln_sfco2inc ) THEN #if defined key_medusa CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) ) CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) ) #elif defined key_hadocc CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) ) #elif defined key_fabm CALL iom_rstput( kt, kt, inum, 'fco2_o3c', pco2_balinc(:,:,:,jp_fabm_o3c) ) CALL iom_rstput( kt, kt, inum, 'fco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) ) #endif ENDIF IF ( ln_pphinc ) THEN #if defined key_medusa CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) ) CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) ) #elif defined key_hadocc CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) ) #elif defined key_fabm CALL iom_rstput( kt, kt, inum, 'ph_o3c', ph_balinc(:,:,:,jp_fabm_o3c) ) CALL iom_rstput( kt, kt, inum, 'ph_o3ba', ph_balinc(:,:,:,jp_fabm_o3ba) ) #endif ENDIF CALL iom_close( inum ) ELSE CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', & & ' Therefore not writing out balancing increments at this timestep', & & ' Something has probably gone wrong somewhere' ) IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// & & TRIM( c_asmbal ) // ' at timestep = ', kt ENDIF END SUBROUTINE asm_bgc_bal_wri !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_bkg_alloc !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_bkg_alloc *** !! !! ** Purpose : allocate time-average arrays for background !! !! ** Method : allocate time-average arrays for background !! !! ** Action : allocate time-average arrays for background !! !! References : asm_bkg_wri !!------------------------------------------------------------------------ !! INTEGER :: ierror !! !!------------------------------------------------------------------------ ALLOCATE( pgrow_avg_tavg(jpi,jpj), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate pgrow_avg_tavg' ) ENDIF pgrow_avg_tavg(:,:) = 0.0 ALLOCATE( ploss_avg_tavg(jpi,jpj), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate ploss_avg_tavg' ) ENDIF ploss_avg_tavg(:,:) = 0.0 ALLOCATE( phyt_avg_tavg(jpi,jpj), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate phyt_avg_tavg' ) ENDIF phyt_avg_tavg(:,:) = 0.0 ALLOCATE( trn_tavg(jpi,jpj,jpk,jptra), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate trn_tavg' ) ENDIF trn_tavg(:,:,:,:) = 0.0 #if defined key_hadocc ALLOCATE( HADOCC_CHL_tavg(jpi,jpj,jpk), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate HADOCC_CHL_tavg' ) ENDIF HADOCC_CHL_tavg(:,:,:) = 0.0 ALLOCATE( cchl_p_tavg(jpi,jpj,jpk), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate cchl_p_tavg' ) ENDIF cchl_p_tavg(:,:,:) = 0.0 #elif defined key_fabm ALLOCATE( totalk_tavg(jpi,jpj,jpk), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate totalk_tavg' ) ENDIF totalk_tavg(:,:,:) = 0.0 #endif END SUBROUTINE asm_bgc_bkg_alloc !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_bkg_tavg(kt, pnumtimes_tavg) !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_bkg_wri *** !! !! ** Purpose : write out bgc background !! !! ** Method : write out bgc background !! !! ** Action : write out bgc background !! !! References : asm_bkg_wri !!------------------------------------------------------------------------ !! INTEGER, INTENT(in ) :: kt ! Current time-step !! REAL(wp), INTENT(in ) :: pnumtimes_tavg ! No of times to average over !! !!------------------------------------------------------------------------ IF (kt == nittrc000) THEN pgrow_avg_tavg(:,:) = 0.0 ploss_avg_tavg(:,:) = 0.0 phyt_avg_tavg(:,:) = 0.0 trn_tavg(:,:,:,:) = 0.0 #if defined key_hadocc HADOCC_CHL_tavg(:,:,:) = 0.0 cchl_p_tavg(:,:,:) = 0.0 #elif defined key_fabm totalk_tavg(:,:,:) = 0.0 #endif ENDIF pgrow_avg_tavg(:,:) = pgrow_avg_tavg(:,:) + pgrow_avg(:,:) / pnumtimes_tavg ploss_avg_tavg(:,:) = ploss_avg_tavg(:,:) + ploss_avg(:,:) / pnumtimes_tavg phyt_avg_tavg(:,:) = phyt_avg_tavg(:,:) + phyt_avg(:,:) / pnumtimes_tavg trn_tavg(:,:,:,:) = trn_tavg(:,:,:,:) + trn(:,:,:,:) / pnumtimes_tavg #if defined key_hadocc HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL_tavg(:,:,:) + HADOCC_CHL(:,:,:) / pnumtimes_tavg cchl_p_tavg(:,:,:) = cchl_p_tavg(:,:,:) + cchl_p(:,:,:) / pnumtimes_tavg #elif defined key_fabm totalk_tavg(:,:,:) = totalk_tavg(:,:,:) + & & fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:) #endif END SUBROUTINE asm_bgc_bkg_tavg !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_bkg_wri( kt, knum, ld_avgbkg ) !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_bkg_wri *** !! !! ** Purpose : write out bgc background !! !! ** Method : write out bgc background !! !! ** Action : write out bgc background !! !! References : asm_bkg_wri !!------------------------------------------------------------------------ !! INTEGER, INTENT(in ) :: kt ! Current time-step INTEGER, INTENT(in ) :: knum ! i/o unit of increments file LOGICAL, INTENT(in ) :: ld_avgbkg ! Averaged background? !! INTEGER :: nitbgcbkg_r ! Period referenced to nit000 !! !!------------------------------------------------------------------------ IF (ld_avgbkg) THEN nitbgcbkg_r = nitavgbkg_r ELSE nitbgcbkg_r = nitbkg_r pgrow_avg_tavg(:,:) = pgrow_avg(:,:) ploss_avg_tavg(:,:) = ploss_avg(:,:) phyt_avg_tavg(:,:) = phyt_avg(:,:) trn_tavg(:,:,:,:) = trn(:,:,:,:) #if defined key_hadocc HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL(:,:,:) cchl_p_tavg(:,:,:) = cchl_p(:,:,:) #elif defined key_fabm totalk_tavg(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:) #endif ENDIF #if defined key_hadocc CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_nut' , trn_tavg(:,:,:,jp_had_nut) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_phy' , trn_tavg(:,:,:,jp_had_phy) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_zoo' , trn_tavg(:,:,:,jp_had_zoo) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_pdn' , trn_tavg(:,:,:,jp_had_pdn) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_dic' , trn_tavg(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_alk' , trn_tavg(:,:,:,jp_had_alk) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_chl' , HADOCC_CHL(:,:,:) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:) ) #elif defined key_medusa CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chn' , trn_tavg(:,:,:,jpchn) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chd' , trn_tavg(:,:,:,jpchd) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phn' , trn_tavg(:,:,:,jpphn) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phd' , trn_tavg(:,:,:,jpphd) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_pds' , trn_tavg(:,:,:,jppds) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zmi' , trn_tavg(:,:,:,jpzmi) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zme' , trn_tavg(:,:,:,jpzme) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_din' , trn_tavg(:,:,:,jpdin) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_sil' , trn_tavg(:,:,:,jpsil) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_fer' , trn_tavg(:,:,:,jpfer) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_det' , trn_tavg(:,:,:,jpdet) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dtc' , trn_tavg(:,:,:,jpdtc) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dic' , trn_tavg(:,:,:,jpdic) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_alk' , trn_tavg(:,:,:,jpalk) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_oxy' , trn_tavg(:,:,:,jpoxy) ) #elif defined key_fabm CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl1' , trn_tavg(:,:,:,jp_fabm_chl1) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl2' , trn_tavg(:,:,:,jp_fabm_chl2) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl3' , trn_tavg(:,:,:,jp_fabm_chl3) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl4' , trn_tavg(:,:,:,jp_fabm_chl4) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1c' , trn_tavg(:,:,:,jp_fabm_p1c) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1n' , trn_tavg(:,:,:,jp_fabm_p1n) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1p' , trn_tavg(:,:,:,jp_fabm_p1p) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1s' , trn_tavg(:,:,:,jp_fabm_p1s) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2c' , trn_tavg(:,:,:,jp_fabm_p2c) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2n' , trn_tavg(:,:,:,jp_fabm_p2n) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2p' , trn_tavg(:,:,:,jp_fabm_p2p) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3c' , trn_tavg(:,:,:,jp_fabm_p3c) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3n' , trn_tavg(:,:,:,jp_fabm_p3n) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3p' , trn_tavg(:,:,:,jp_fabm_p3p) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4c' , trn_tavg(:,:,:,jp_fabm_p4c) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4n' , trn_tavg(:,:,:,jp_fabm_p4n) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4p' , trn_tavg(:,:,:,jp_fabm_p4p) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z4c' , trn_tavg(:,:,:,jp_fabm_z4c) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5c' , trn_tavg(:,:,:,jp_fabm_z5c) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5n' , trn_tavg(:,:,:,jp_fabm_z5n) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5p' , trn_tavg(:,:,:,jp_fabm_z5p) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6c' , trn_tavg(:,:,:,jp_fabm_z6c) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6n' , trn_tavg(:,:,:,jp_fabm_z6n) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6p' , trn_tavg(:,:,:,jp_fabm_z6p) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n1p' , trn_tavg(:,:,:,jp_fabm_n1p) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n3n' , trn_tavg(:,:,:,jp_fabm_n3n) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n4n' , trn_tavg(:,:,:,jp_fabm_n4n) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n5s' , trn_tavg(:,:,:,jp_fabm_n5s) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o2o' , trn_tavg(:,:,:,jp_fabm_o2o) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3c' , trn_tavg(:,:,:,jp_fabm_o3c) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ba' , trn_tavg(:,:,:,jp_fabm_o3ba) ) CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ta' , totalk_tavg ) #endif END SUBROUTINE asm_bgc_bkg_wri !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog ) !!------------------------------------------------------------------------ !! *** ROUTINE asm_bgc_init_incs *** !! !! ** Purpose : convert log increments to non-log !! !! ** Method : need to account for model background, !! cannot simply do 10^log_inc. Need to: !! 1) Add log_inc to log10(background) to get log10(analysis) !! 2) Take 10^log10(analysis) to get analysis !! 3) Subtract background from analysis to get non-log incs !! !! ** Action : return non-log increments !! !! References : !!------------------------------------------------------------------------ !! REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pbkg ! Background REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pinc_log ! Log incs REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs ! INTEGER :: ji, jj ! Loop counters !! !!------------------------------------------------------------------------ DO jj = 1, jpj DO ji = 1, jpi IF ( pbkg(ji,jj) > 0.0 ) THEN pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + & & pinc_log(ji,jj) ) & & - pbkg(ji,jj) ELSE pinc_nonlog(ji,jj) = 0.0 ENDIF END DO END DO END SUBROUTINE asm_bgc_unlog_2d !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau ) !!------------------------------------------------------------------------ !! *** ROUTINE phyto2d_asm_inc *** !! !! ** Purpose : Apply the chlorophyll assimilation increments. !! !! ** Method : Calculate increments to state variables using nitrogen !! balancing. !! Direct initialization or Incremental Analysis Updating. !! !! ** Action : !!------------------------------------------------------------------------ INTEGER, INTENT(IN) :: kt ! Current time step LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights ! INTEGER :: it ! Index REAL(wp) :: zincwgt ! IAU weight for current time step REAL(wp) :: zincper ! IAU interval in seconds REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg #if defined key_medusa REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg #elif defined key_fabm REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldin ! Local chldin incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldin ! Local chldin bkg REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnan ! Local chlnan incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnan ! Local chlnan bkg REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlpic ! Local chlpic incs REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlpic ! Local chlpic bkg #endif !!------------------------------------------------------------------------ IF ( kt <= nit000 ) THEN ! Un-log any log increments for passing to balancing routines ! Remember that two sets of non-log increments should not be ! expected to be in the same ratio as their log equivalents ! Total chlorophyll IF ( ln_slchltotinc ) THEN #if defined key_medusa zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd) #elif defined key_hadocc zbkg_chltot(:,:) = chl_bkg(:,:,1) #elif defined key_fabm zbkg_chltot(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1) + & & tracer_bkg(:,:,1,jp_fabm_chl2) + & & tracer_bkg(:,:,1,jp_fabm_chl3) + & & tracer_bkg(:,:,1,jp_fabm_chl4) #endif CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot ) ELSE IF ( ln_schltotinc ) THEN zinc_chltot(:,:) = schltot_bkginc(:,:) ELSE zinc_chltot(:,:) = 0.0 ENDIF #if defined key_medusa || defined key_fabm ! Diatom chlorophyll IF ( ln_slchldiainc ) THEN #if defined key_medusa zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd) #elif defined key_fabm zbkg_chldia(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1) #endif CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia ) ELSE zinc_chldia(:,:) = 0.0 ENDIF #endif #if defined key_medusa ! Non-diatom chlorophyll IF ( ln_slchlnoninc ) THEN zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn) CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon ) ELSE zinc_chlnon(:,:) = 0.0 ENDIF #endif #if defined key_fabm ! Nanophytoplankton chlorophyll IF ( ln_slchlnaninc ) THEN zbkg_chlnan(:,:) = tracer_bkg(:,:,1,jp_fabm_chl2) CALL asm_bgc_unlog_2d( zbkg_chlnan, slchlnan_bkginc, zinc_chlnan ) ELSE zinc_chlnan(:,:) = 0.0 ENDIF ! Picophytoplankton chlorophyll IF ( ln_slchlpicinc ) THEN zbkg_chlpic(:,:) = tracer_bkg(:,:,1,jp_fabm_chl3) CALL asm_bgc_unlog_2d( zbkg_chlpic, slchlpic_bkginc, zinc_chlpic ) ELSE zinc_chlpic(:,:) = 0.0 ENDIF ! Dinoflagellate chlorophyll IF ( ln_slchldininc ) THEN zbkg_chldin(:,:) = tracer_bkg(:,:,1,jp_fabm_chl4) CALL asm_bgc_unlog_2d( zbkg_chldin, slchldin_bkginc, zinc_chldin ) ELSE zinc_chldin(:,:) = 0.0 ENDIF #endif ! Total phytoplankton carbon IF ( ln_slphytotinc ) THEN #if defined key_medusa zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) #elif defined key_hadocc zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p #endif CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot ) ELSE zinc_phytot(:,:) = 0.0 ENDIF #if defined key_medusa ! Diatom phytoplankton carbon IF ( ln_slphydiainc ) THEN zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia ) ELSE zinc_phydia(:,:) = 0.0 ENDIF #endif #if defined key_medusa ! Non-diatom phytoplankton carbon IF ( ln_slphynoninc ) THEN zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon ) ELSE zinc_phynon(:,:) = 0.0 ENDIF #endif ! Select mixed layer IF ( ll_asmdin ) THEN #if defined key_top && ( defined key_hadocc || defined key_medusa || defined key_fabm ) CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', & & ' Mixed layer depth taken to be background maximum mld_max_bkg' ) zmld(:,:) = mld_max_bkg(:,:) #else CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' ) #endif ELSE SELECT CASE( mld_choice_bgc ) CASE ( 1 ) ! Turbocline/mixing depth [W points] zmld(:,:) = hmld(:,:) CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points] zmld(:,:) = hmlp(:,:) CASE ( 3 ) ! Kara MLD [Interpolated] #if defined key_karaml IF ( ln_kara ) THEN zmld(:,:) = hmld_kara(:,:) ELSE CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', & & ' but ln_kara=.false.' ) ENDIF #else CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', & & ' but is not defined' ) #endif CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points] zmld(:,:) = hmld_tref(:,:) CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points] zmld(:,:) = hmlpt(:,:) END SELECT ENDIF zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt #if defined key_medusa CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), & & zinc_chltot, & & ln_slchldiainc, & & zinc_chldia, & & ln_slchlnoninc, & & zinc_chlnon, & & ln_slphytotinc, & & zinc_phytot, & & ln_slphydiainc, & & zinc_phydia, & & ln_slphynoninc, & & zinc_phynon, & & zincper, & & rn_maxchlinc, ln_phytobal, zmld, & & pgrow_avg_bkg, ploss_avg_bkg, & & phyt_avg_bkg, mld_max_bkg, & & tracer_bkg, phyto2d_balinc ) #elif defined key_hadocc CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), & & zinc_chltot, & & ln_slphytotinc, & & zinc_phytot, & & zincper, & & rn_maxchlinc, ln_phytobal, zmld, & & pgrow_avg_bkg, ploss_avg_bkg, & & phyt_avg_bkg, mld_max_bkg, & & cchl_p_bkg(:,:,1), & & tracer_bkg, phyto2d_balinc ) #elif defined key_fabm CALL asm_phyto2d_bal_ersem( (ln_slchltotinc .OR. ln_schltotinc), & & zinc_chltot, & & ln_slchldiainc, & & zinc_chldia, & & ln_slchlnaninc, & & zinc_chlnan, & & ln_slchlpicinc, & & zinc_chlpic, & & ln_slchldininc, & & zinc_chldin, & & zincper, & & rn_maxchlinc, ln_phytobal, zmld, & & pgrow_avg_bkg, ploss_avg_bkg, & & phyt_avg_bkg, mld_max_bkg, & & totalk_bkg, & & tracer_bkg, phyto2d_balinc ) #else CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', & & 'but not defined a biogeochemical model' ) #endif ENDIF IF ( ll_asmiau ) THEN !-------------------------------------------------------------------- ! Incremental Analysis Updating !-------------------------------------------------------------------- IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN it = kt - nit000 + 1 zincwgt = pwgtiau(it) ! IAU weight for the current time step ! note this is not a tendency so should not be divided by rdt IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', & & kt,' with IAU weight = ', pwgtiau(it) WRITE(numout,*) '~~~~~~~~~~~~' ENDIF ! Update the biogeochemical variables ! Add directly to trn and trb, rather than to tra, because tra gets ! reset to zero at the start of trc_stp, called after this routine #if defined key_medusa WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. & & trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + & & phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt END WHERE #elif defined key_hadocc WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. & & trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & & phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + & & phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt END WHERE #elif defined key_fabm WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & & phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + & & phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt END WHERE #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ELSEIF ( ll_asmdin ) THEN !-------------------------------------------------------------------- ! Direct Initialization !-------------------------------------------------------------------- IF ( kt == nitdin_r ) THEN neuler = 0 ! Force Euler forward step ! Initialize the now fields with the background + increment ! Background currently is what the model is initialised with CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', & & ' Background state is taken from model rather than background file' ) #if defined key_medusa WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. & & trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) END WHERE #elif defined key_hadocc WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. & & trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp ) trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & & phyto2d_balinc(:,:,:,jp_had0:jp_had1) trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) END WHERE #elif defined key_fabm WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp ) trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & & phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) END WHERE #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF ! END SUBROUTINE phyto2d_asm_inc !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau ) !!------------------------------------------------------------------------ !! *** ROUTINE phyto3d_asm_inc *** !! !! ** Purpose : Apply the profile chlorophyll assimilation increments. !! !! ** Method : Calculate increments to state variables. !! Direct initialization or Incremental Analysis Updating. !! !! ** Action : !!------------------------------------------------------------------------ INTEGER, INTENT(IN) :: kt ! Current time step LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights ! INTEGER :: ji, jj, jk ! Loop counters INTEGER :: it ! Index REAL(wp) :: zincwgt ! IAU weight for timestep REAL(wp) :: zfrac_chn ! Fraction of jpchn REAL(wp) :: zfrac_chd ! Fraction of jpchd REAL(wp) :: zfrac_chl1 ! Fraction of jp_fabm_chl1 REAL(wp) :: zfrac_chl2 ! Fraction of jp_fabm_chl2 REAL(wp) :: zfrac_chl3 ! Fraction of jp_fabm_chl3 REAL(wp) :: zfrac_chl4 ! Fraction of jp_fabm_chl4 REAL(wp) :: zrat_phn_chn ! jpphn:jpchn ratio REAL(wp) :: zrat_phd_chd ! jpphd:jpchd ratio REAL(wp) :: zrat_pds_chd ! jppds:jpchd ratio REAL(wp) :: zrat_p1c_chl1 ! jp_fabm_p1c:jp_fabm_chl1 ratio REAL(wp) :: zrat_p1n_chl1 ! jp_fabm_p1n:jp_fabm_chl1 ratio REAL(wp) :: zrat_p1p_chl1 ! jp_fabm_p1p:jp_fabm_chl1 ratio REAL(wp) :: zrat_p1s_chl1 ! jp_fabm_p1s:jp_fabm_chl1 ratio REAL(wp) :: zrat_p2c_chl2 ! jp_fabm_p2c:jp_fabm_chl2 ratio REAL(wp) :: zrat_p2n_chl2 ! jp_fabm_p2n:jp_fabm_chl2 ratio REAL(wp) :: zrat_p2p_chl2 ! jp_fabm_p2p:jp_fabm_chl2 ratio REAL(wp) :: zrat_p3c_chl3 ! jp_fabm_p3c:jp_fabm_chl3 ratio REAL(wp) :: zrat_p3n_chl3 ! jp_fabm_p3n:jp_fabm_chl3 ratio REAL(wp) :: zrat_p3p_chl3 ! jp_fabm_p3p:jp_fabm_chl3 ratio REAL(wp) :: zrat_p4c_chl4 ! jp_fabm_p4c:jp_fabm_chl4 ratio REAL(wp) :: zrat_p4n_chl4 ! jp_fabm_p4n:jp_fabm_chl4 ratio REAL(wp) :: zrat_p4p_chl4 ! jp_fabm_p4p:jp_fabm_chl4 ratio REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc ! Chlorophyll increments REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl ! Chlorophyll background !!------------------------------------------------------------------------ IF ( kt <= nit000 ) THEN IF ( ln_plchltotinc ) THEN ! Convert log10(chlorophyll) increment back to a chlorophyll increment ! In order to transform logchl incs to chl incs, need to account for model ! background, cannot simply do 10^logchl_bkginc. Need to: ! 1) Add logchl inc to log10(background) to get log10(analysis) ! 2) Take 10^log10(analysis) to get analysis ! 3) Subtract background from analysis to get chl incs ! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value #if defined key_medusa bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd) #elif defined key_hadocc bkg_chl(:,:,:) = chl_bkg(:,:,:) #elif defined key_fabm bkg_chl(:,:,:) = tracer_bkg(:,:,:,jp_fabm_chl1) + & & tracer_bkg(:,:,:,jp_fabm_chl2) + & & tracer_bkg(:,:,:,jp_fabm_chl3) + & & tracer_bkg(:,:,:,jp_fabm_chl4) #endif DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk) IF ( rn_maxchlinc > 0.0 ) THEN chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) ) ENDIF ELSE chl_inc(ji,jj,jk) = 0.0 ENDIF END DO END DO END DO ELSE IF ( ln_pchltotinc ) THEN DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi IF ( rn_maxchlinc > 0.0 ) THEN chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) ) ELSE chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk) ENDIF END DO END DO END DO ENDIF #if defined key_medusa ! Loop over each grid point partioning the increments based on existing ratios DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd)) zfrac_chd = 1.0 - zfrac_chn phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn) zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd) zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd) phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd ENDIF END DO END DO END DO #elif defined key_hadocc phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:) #elif defined key_fabm ! Loop over each grid point partioning the increments based on existing ratios DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_chl1) > 0.0 ) .AND. & & ( tracer_bkg(ji,jj,jk,jp_fabm_chl2) > 0.0 ) .AND. & & ( tracer_bkg(ji,jj,jk,jp_fabm_chl3) > 0.0 ) .AND. & & ( tracer_bkg(ji,jj,jk,jp_fabm_chl4) > 0.0 ) ) THEN zfrac_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_chl1) / bkg_chl(ji,jj,jk) zfrac_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_chl2) / bkg_chl(ji,jj,jk) zfrac_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_chl3) / bkg_chl(ji,jj,jk) zfrac_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_chl4) / bkg_chl(ji,jj,jk) phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) = chl_inc(ji,jj,jk) * zfrac_chl1 phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) = chl_inc(ji,jj,jk) * zfrac_chl2 phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) = chl_inc(ji,jj,jk) * zfrac_chl3 phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) = chl_inc(ji,jj,jk) * zfrac_chl4 zrat_p1c_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1c) / tracer_bkg(ji,jj,jk,jp_fabm_chl1) zrat_p1n_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1n) / tracer_bkg(ji,jj,jk,jp_fabm_chl1) zrat_p1p_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1p) / tracer_bkg(ji,jj,jk,jp_fabm_chl1) zrat_p1s_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1s) / tracer_bkg(ji,jj,jk,jp_fabm_chl1) zrat_p2c_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2c) / tracer_bkg(ji,jj,jk,jp_fabm_chl2) zrat_p2n_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2n) / tracer_bkg(ji,jj,jk,jp_fabm_chl2) zrat_p2p_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2p) / tracer_bkg(ji,jj,jk,jp_fabm_chl2) zrat_p3c_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3c) / tracer_bkg(ji,jj,jk,jp_fabm_chl3) zrat_p3n_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3n) / tracer_bkg(ji,jj,jk,jp_fabm_chl3) zrat_p3p_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3p) / tracer_bkg(ji,jj,jk,jp_fabm_chl3) zrat_p4c_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4c) / tracer_bkg(ji,jj,jk,jp_fabm_chl4) zrat_p4n_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4n) / tracer_bkg(ji,jj,jk,jp_fabm_chl4) zrat_p4p_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4p) / tracer_bkg(ji,jj,jk,jp_fabm_chl4) phyto3d_balinc(ji,jj,jk,jp_fabm_p1c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1c_chl1 phyto3d_balinc(ji,jj,jk,jp_fabm_p1n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1n_chl1 phyto3d_balinc(ji,jj,jk,jp_fabm_p1p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1p_chl1 phyto3d_balinc(ji,jj,jk,jp_fabm_p1s) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1s_chl1 phyto3d_balinc(ji,jj,jk,jp_fabm_p2c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2c_chl2 phyto3d_balinc(ji,jj,jk,jp_fabm_p2n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2n_chl2 phyto3d_balinc(ji,jj,jk,jp_fabm_p2p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2p_chl2 phyto3d_balinc(ji,jj,jk,jp_fabm_p3c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3c_chl3 phyto3d_balinc(ji,jj,jk,jp_fabm_p3n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3n_chl3 phyto3d_balinc(ji,jj,jk,jp_fabm_p3p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3p_chl3 phyto3d_balinc(ji,jj,jk,jp_fabm_p4c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4c_chl4 phyto3d_balinc(ji,jj,jk,jp_fabm_p4n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4n_chl4 phyto3d_balinc(ji,jj,jk,jp_fabm_p4p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4p_chl4 ENDIF END DO END DO END DO #else CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', & & 'but not defined a biogeochemical model' ) #endif ENDIF IF ( ll_asmiau ) THEN !-------------------------------------------------------------------- ! Incremental Analysis Updating !-------------------------------------------------------------------- IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN it = kt - nit000 + 1 zincwgt = pwgtiau(it) ! IAU weight for the current time step ! note this is not a tendency so should not be divided by rdt IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', & & kt,' with IAU weight = ', pwgtiau(it) WRITE(numout,*) '~~~~~~~~~~~~' ENDIF ! Update the biogeochemical variables ! Add directly to trn and trb, rather than to tra, because tra gets ! reset to zero at the start of trc_stp, called after this routine #if defined key_medusa WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. & & trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + & & phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt END WHERE #elif defined key_hadocc WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. & & trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & & phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + & & phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt END WHERE #elif defined key_fabm WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & & phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + & & phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt END WHERE #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ELSEIF ( ll_asmdin ) THEN !-------------------------------------------------------------------- ! Direct Initialization !-------------------------------------------------------------------- IF ( kt == nitdin_r ) THEN neuler = 0 ! Force Euler forward step ! Initialize the now fields with the background + increment ! Background currently is what the model is initialised with CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', & & ' Background state is taken from model rather than background file' ) #if defined key_medusa WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. & & trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) END WHERE #elif defined key_hadocc WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. & & trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp ) trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & & phyto3d_balinc(:,:,:,jp_had0:jp_had1) trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) END WHERE #elif defined key_fabm WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp ) trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & & phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) END WHERE #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF ! END SUBROUTINE phyto3d_asm_inc !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, & & ll_trainc, pt_bkginc, ps_bkginc ) !!------------------------------------------------------------------------ !! *** ROUTINE pco2_asm_inc *** !! !! ** Purpose : Apply the pco2/fco2 assimilation increments. !! !! ** Method : Calculate increments to state variables using carbon !! balancing. !! Direct initialization or Incremental Analysis Updating. !! !! ** Action : !!------------------------------------------------------------------------ INTEGER, INTENT(IN) :: kt ! Current time step LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments ! INTEGER :: ji, jj, jk ! Loop counters INTEGER :: it ! Index INTEGER :: jkmax ! Index of mixed layer depth REAL(wp) :: zincwgt ! IAU weight for current time step REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure ! Coefficients for fCO2 to pCO2 conversion REAL(wp) :: zcoef_fco2_1 = -1636.75 REAL(wp) :: zcoef_fco2_2 = 12.0408 REAL(wp) :: zcoef_fco2_3 = 0.0327957 REAL(wp) :: zcoef_fco2_4 = 0.0000316528 REAL(wp) :: zcoef_fco2_5 = 2.0 REAL(wp) :: zcoef_fco2_6 = 57.7 REAL(wp) :: zcoef_fco2_7 = 0.118 REAL(wp) :: zcoef_fco2_8 = 82.0578 !!------------------------------------------------------------------------ IF ( kt <= nit000 ) THEN !---------------------------------------------------------------------- ! DIC and alkalinity balancing !---------------------------------------------------------------------- IF ( ln_sfco2inc ) THEN #if defined key_medusa && defined key_roam ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine patm(1) = 1.0 phyd(1) = 0.0 DO jj = 1, jpj DO ji = 1, jpi CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) ) END DO END DO #else ! If assimilating fCO2, then convert to pCO2 using temperature ! See flux_gas.F90 within HadOCC for details of calculation pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) / & & EXP((zcoef_fco2_1 + & & zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - & & zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + & & zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + & & zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / & & (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0))) #endif ELSE pco2_bkginc_temp(:,:) = spco2_bkginc(:,:) ENDIF ! Account for physics assimilation if required IF ( ll_trainc ) THEN tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1) sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1) ELSE tem_bkg_temp(:,:) = tsn(:,:,1,1) sal_bkg_temp(:,:) = tsn(:,:,1,2) ENDIF #if defined key_medusa ! Account for phytoplankton balancing if required IF ( ln_phytobal ) THEN dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic) alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk) ELSE dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) ENDIF CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & & tem_bkg_temp(:,:), sal_bkg_temp(:,:), & & pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) ) #elif defined key_hadocc ! Account for phytoplankton balancing if required IF ( ln_phytobal ) THEN dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic) alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk) ELSE dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) ENDIF CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & & tem_bkg_temp(:,:), sal_bkg_temp(:,:), & & pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) ) #elif defined key_fabm ! Account for phytoplankton balancing if required IF ( ln_phytobal ) THEN dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c) + phyto2d_balinc(:,:,1,jp_fabm_o3c) alk_bkg_temp(:,:) = totalk_bkg(:,:,1) + phyto2d_balinc(:,:,1,jp_fabm_o3ba) ELSE dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c) alk_bkg_temp(:,:) = totalk_bkg(:,:,1) ENDIF CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & & tem_bkg_temp(:,:), sal_bkg_temp(:,:), & & pco2_balinc(:,:,1,jp_fabm_o3c), pco2_balinc(:,:,1,jp_fabm_o3ba) ) #else CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', & & 'but not defined a biogeochemical model' ) #endif ! Select mixed layer IF ( ll_asmdin ) THEN #if defined key_hadocc || defined key_medusa || defined key_fabm CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', & & ' Mixed layer depth taken to be background maximum mld_max_bkg' ) zmld(:,:) = mld_max_bkg(:,:) #else CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' ) #endif ELSE SELECT CASE( mld_choice_bgc ) CASE ( 1 ) ! Turbocline/mixing depth [W points] zmld(:,:) = hmld(:,:) CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points] zmld(:,:) = hmlp(:,:) CASE ( 3 ) ! Kara MLD [Interpolated] #if defined key_karaml IF ( ln_kara ) THEN zmld(:,:) = hmld_kara(:,:) ELSE CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & & ' but ln_kara=.false.' ) ENDIF #else CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & & ' but is not defined' ) #endif CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points] zmld(:,:) = hmld_tref(:,:) CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points] zmld(:,:) = hmlpt(:,:) END SELECT ENDIF ! Propagate through mixed layer DO jj = 1, jpj DO ji = 1, jpi ! jkmax = jpk-1 DO jk = jpk-1, 1, -1 IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. & & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN zmld(ji,jj) = gdepw_n(ji,jj,jk+1) jkmax = jk ENDIF END DO ! #if defined key_top DO jk = 2, jkmax pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:) END DO #endif ! END DO END DO ENDIF IF ( ll_asmiau ) THEN !-------------------------------------------------------------------- ! Incremental Analysis Updating !-------------------------------------------------------------------- IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN it = kt - nit000 + 1 zincwgt = pwgtiau(it) ! IAU weight for the current time step ! note this is not a tendency so should not be divided by rdt IF(lwp) THEN WRITE(numout,*) IF ( ln_spco2inc ) THEN WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', & & kt,' with IAU weight = ', pwgtiau(it) ELSE IF ( ln_sfco2inc ) THEN WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', & & kt,' with IAU weight = ', pwgtiau(it) ENDIF WRITE(numout,*) '~~~~~~~~~~~~' ENDIF ! Update the biogeochemical variables ! Add directly to trn and trb, rather than to tra, because tra gets ! reset to zero at the start of trc_stp, called after this routine #if defined key_medusa WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. & & trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + & & pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt END WHERE #elif defined key_hadocc WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. & & trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & & pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + & & pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt END WHERE #elif defined key_fabm WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & & pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + & & pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt END WHERE #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ELSEIF ( ll_asmdin ) THEN !-------------------------------------------------------------------- ! Direct Initialization !-------------------------------------------------------------------- IF ( kt == nitdin_r ) THEN neuler = 0 ! Force Euler forward step ! Initialize the now fields with the background + increment ! Background currently is what the model is initialised with CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', & & ' Background state is taken from model rather than background file' ) #if defined key_medusa WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. & & trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & pco2_balinc(:,:,:,jp_msa0:jp_msa1) trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) END WHERE #elif defined key_hadocc WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. & & trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp ) trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & & pco2_balinc(:,:,:,jp_had0:jp_had1) trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) END WHERE #elif defined key_fabm WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp ) trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & & pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) END WHERE #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF ! END SUBROUTINE pco2_asm_inc !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, & & ll_trainc, pt_bkginc, ps_bkginc ) !!------------------------------------------------------------------------ !! *** ROUTINE ph_asm_inc *** !! !! ** Purpose : Apply the pH assimilation increments. !! !! ** Method : Calculate increments to DIC and alkalinity from pH !! Use a similar approach to the pCO2 scheme !! Direct initialization or Incremental Analysis Updating. !! !! ** Action : !!------------------------------------------------------------------------ INTEGER, INTENT(IN) :: kt ! Current time step LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments REAL(wp) :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation REAL(wp) :: DIC_IN, ALK_IN ! DIC/alk in pH calculation REAL(wp) :: PH_OUT ! pH from pH calculation REAL(wp) :: ph_bkg ! pH from pH calculation REAL(wp) :: ph_dic ! pH from pH calculation REAL(wp) :: ph_alk ! pH from pH calculation REAL(wp) :: dph_ddic ! Change in pH wrt DIC REAL(wp) :: dph_dalk ! Change in pH wrt alk REAL(wp) :: weight ! Increment weighting REAL(wp) :: zincwgt ! IAU weight for current time step REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp ! DIC background REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp ! Alkalinity background REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp ! DIN background REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp ! Silicate background REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp ! Temperature background REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp ! Salinity background REAL(wp), DIMENSION(19) :: dummy_out ! Dummy variables INTEGER :: ji, jj, jk, jx ! Loop counters INTEGER :: it ! Index !!------------------------------------------------------------------------ #if ! defined key_medusa CALL ctl_stop( ' pH balancing only implemented for MEDUSA' ) #else IF ( kt <= nit000 ) THEN !---------------------------------------------------------------------- ! DIC and alkalinity balancing !---------------------------------------------------------------------- ! Account for physics assimilation if required IF ( ll_trainc ) THEN tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:) sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:) ELSE tem_bkg_temp(:,:,:) = tsn(:,:,:,1) sal_bkg_temp(:,:,:) = tsn(:,:,:,2) ENDIF ! Account for phytoplankton balancing if required IF ( ln_phytobal ) THEN dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic) alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk) din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin) sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil) ELSE dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) ENDIF ! Account for pCO2 balancing if required IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic) alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk) ENDIF ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk ! This requires three calls to the MOCSY carbonate package ! One to find pH at background DIC and alkalinity ! One to find pH when DIC is increased by 10 ! One to find pH when alkalinity is increased by 10 ! Then calculate the gradients DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN DO jx = 1, 3 IF ( jx == 1 ) THEN DIC_IN = dic_bkg_temp(ji,jj,jk) ALK_IN = alk_bkg_temp(ji,jj,jk) ELSE IF ( jx == 2 ) THEN DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch ALK_IN = alk_bkg_temp(ji,jj,jk) ELSE IF ( jx == 3 ) THEN DIC_IN = dic_bkg_temp(ji,jj,jk) ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch ENDIF CALL mocsy_interface(tsn(ji,jj,jk,jp_tem), & ! IN: temperature & tsn(ji,jj,jk,jp_sal), & ! salinity & ALK_IN, & ! alkalinity & DIC_IN, & ! DIC & sil_bkg_temp(ji,jj,jk), & ! silicate & din_bkg_temp(ji,jj,jk)/16.0, & ! phosphate (Redfield from DIN) & 1.0, & ! atmospheric pressure (dummy) & fsdept(ji,jj,jk), & ! depth & gphit(ji,jj), & ! latitude & 1.0, & ! gas transfer (dummy) & 1.0, & ! atmospheric xCO2 (dummy) & 1, & ! number of points & PH_OUT, & ! OUT: pH & dummy_out(1), dummy_out(2), & ! stuff we don't care about & dummy_out(3), dummy_out(4), & & dummy_out(5), dummy_out(6), & & dummy_out(7), dummy_out(8), & & dummy_out(9), dummy_out(10), & & dummy_out(11), dummy_out(12), & & dummy_out(13), dummy_out(14), & & dummy_out(15), dummy_out(16), & & dummy_out(17), dummy_out(18), & & dummy_out(19) ) IF ( jx == 1 ) THEN ph_bkg = PH_OUT ELSE IF ( jx == 2 ) THEN ph_dic = PH_OUT ELSE IF ( jx == 3 ) THEN ph_alk = PH_OUT ENDIF END DO dph_ddic = (ph_dic - ph_bkg) / zsearch dph_dalk = (ph_alk - ph_bkg) / zsearch weight = pph_bkginc(ji,jj,jk) / ( (dph_ddic*dph_ddic) + (dph_dalk*dph_dalk) ) ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk ENDIF END DO END DO END DO ENDIF IF ( ll_asmiau ) THEN !-------------------------------------------------------------------- ! Incremental Analysis Updating !-------------------------------------------------------------------- IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN it = kt - nit000 + 1 zincwgt = pwgtiau(it) ! IAU weight for the current time step ! note this is not a tendency so should not be divided by rdt IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', & & kt,' with IAU weight = ', pwgtiau(it) WRITE(numout,*) '~~~~~~~~~~~~' ENDIF ! Update the biogeochemical variables ! Add directly to trn and trb, rather than to tra, because tra gets ! reset to zero at the start of trc_stp, called after this routine WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. & & trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + & & ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt END WHERE ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ELSEIF ( ll_asmdin ) THEN !-------------------------------------------------------------------- ! Direct Initialization !-------------------------------------------------------------------- IF ( kt == nitdin_r ) THEN neuler = 0 ! Force Euler forward step ! Initialize the now fields with the background + increment ! Background currently is what the model is initialised with CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', & & ' Background state is taken from model rather than background file' ) WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. & & trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & ph_balinc(:,:,:,jp_msa0:jp_msa1) trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) END WHERE ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF #endif ! END SUBROUTINE ph_asm_inc !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_asm_inc *** !! !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments. !! !! ** Method : Direct initialization or Incremental Analysis Updating. !! !! ** Action : !!---------------------------------------------------------------------- INTEGER, INTENT(IN) :: kt ! Current time step LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights ! INTEGER :: it ! Index REAL(wp) :: zincwgt ! IAU weight for current time step !!---------------------------------------------------------------------- IF ( kt <= nit000 ) THEN !---------------------------------------------------------------------- ! Remove any other balancing increments !---------------------------------------------------------------------- IF ( ln_pno3inc ) THEN #if defined key_hadocc || defined key_medusa || defined key_fabm #if defined key_hadocc it = jp_had_nut #elif defined key_medusa it = jpdin #elif defined key_fabm it = jp_fabm_n3n #endif IF ( ln_phytobal ) THEN pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) ENDIF IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it) ENDIF IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it) ENDIF IF ( ln_pphinc ) THEN pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it) ENDIF #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_psi4inc ) THEN #if defined key_medusa || defined key_fabm #if defined key_medusa it = jpsil #elif defined key_fabm it = jp_fabm_n5s #endif IF ( ln_phytobal ) THEN psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) ENDIF IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it) ENDIF IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it) ENDIF IF ( ln_pphinc ) THEN psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it) ENDIF #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_pdicinc ) THEN #if defined key_hadocc || defined key_medusa || defined key_fabm #if defined key_hadocc it = jp_had_dic #elif defined key_medusa it = jpdic #elif defined key_fabm it = jp_fabm_o3c #endif IF ( ln_phytobal ) THEN pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) ENDIF IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it) ENDIF IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it) ENDIF IF ( ln_pphinc ) THEN pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it) ENDIF #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_palkinc ) THEN #if defined key_hadocc || defined key_medusa || defined key_fabm #if defined key_hadocc it = jp_had_alk #elif defined key_medusa it = jpalk #elif defined key_fabm it = jp_fabm_o3ba #endif IF ( ln_phytobal ) THEN palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) ENDIF IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it) ENDIF IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it) ENDIF IF ( ln_pphinc ) THEN palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it) ENDIF #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_po2inc ) THEN #if defined key_medusa || defined key_fabm #if defined key_medusa it = jpoxy #elif defined key_fabm it = jp_fabm_o2o #endif IF ( ln_phytobal ) THEN po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) ENDIF IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it) ENDIF IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it) ENDIF IF ( ln_pphinc ) THEN po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it) ENDIF #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF ENDIF IF ( ll_asmiau ) THEN !-------------------------------------------------------------------- ! Incremental Analysis Updating !-------------------------------------------------------------------- IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN it = kt - nit000 + 1 zincwgt = pwgtiau(it) ! IAU weight for the current time step ! note this is not a tendency so should not be divided by rdt IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', & & kt,' with IAU weight = ', pwgtiau(it) WRITE(numout,*) '~~~~~~~~~~~~' ENDIF ! Update the 3D BGC variables ! Add directly to trn and trb, rather than to tra, because tra gets ! reset to zero at the start of trc_stp, called after this routine ! Don't apply increments if they'll take concentrations negative IF ( ln_pno3inc ) THEN #if defined key_hadocc WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt END WHERE #elif defined key_medusa WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt END WHERE #elif defined key_fabm WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt trb(:,:,:,jp_fabm_n3n) = trb(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_psi4inc ) THEN #if defined key_medusa WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt END WHERE #elif defined key_fabm WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt trb(:,:,:,jp_fabm_n5s) = trb(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_pdicinc ) THEN #if defined key_hadocc WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt END WHERE #elif defined key_medusa WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt END WHERE #elif defined key_fabm WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt trb(:,:,:,jp_fabm_o3c) = trb(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_palkinc ) THEN #if defined key_hadocc WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt END WHERE #elif defined key_medusa WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt END WHERE #elif defined key_fabm WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt trb(:,:,:,jp_fabm_o3ba) = trb(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_po2inc ) THEN #if defined key_medusa WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt END WHERE #elif defined key_fabm WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp ) trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt trb(:,:,:,jp_fabm_o2o) = trb(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( kt == nitiaufin_r ) THEN IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc ) IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc ) IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc ) IF ( ln_palkinc ) DEALLOCATE( palk_bkginc ) IF ( ln_po2inc ) DEALLOCATE( po2_bkginc ) ENDIF ENDIF ELSEIF ( ll_asmdin ) THEN !-------------------------------------------------------------------- ! Direct Initialization !-------------------------------------------------------------------- IF ( kt == nitdin_r ) THEN neuler = 0 ! Force Euler forward step ! Initialize the now fields with the background + increment ! Background currently is what the model is initialised with CALL ctl_warn( ' Doing direct initialisation with 3D BGC assimilation', & & ' Background state is taken from model rather than background file' ) IF ( ln_pno3inc ) THEN #if defined key_hadocc WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) END WHERE #elif defined key_medusa WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) trb(:,:,:,jpdin) = trn(:,:,:,jpdin) END WHERE #elif defined key_fabm WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) trb(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_psi4inc ) THEN #if defined key_medusa WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) trb(:,:,:,jpsil) = trn(:,:,:,jpsil) END WHERE #elif defined key_fabm WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) trb(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_pdicinc ) THEN #if defined key_hadocc WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) END WHERE #elif defined key_medusa WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) trb(:,:,:,jpdic) = trn(:,:,:,jpdic) END WHERE #elif defined key_fabm WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) trb(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_palkinc ) THEN #if defined key_hadocc WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) END WHERE #elif defined key_medusa WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) trb(:,:,:,jpalk) = trn(:,:,:,jpalk) END WHERE #elif defined key_fabm WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) trb(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_po2inc ) THEN #if defined key_medusa WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy) END WHERE #elif defined key_fabm WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. & & trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) > 0.0_wp ) trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) trb(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) END WHERE #else CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) #endif ENDIF IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc ) IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc ) IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc ) IF ( ln_palkinc ) DEALLOCATE( palk_bkginc ) IF ( ln_po2inc ) DEALLOCATE( po2_bkginc ) ENDIF ! ENDIF ! END SUBROUTINE bgc3d_asm_inc !!=========================================================================== END MODULE asmbgc