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_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, & & hmlp, & & hmlpt USE asmpar, ONLY: & ! assimilation parameters & c_asmbkg, & & c_asmbal, & & nitbkg_r, & & nitdin_r, & & nitiaustr_r, & & nitiaufin_r #if defined key_top USE trc, ONLY: & ! passive tracer variables & trn, & & trb USE par_trc, ONLY: & ! passive tracer parameters & jptra #endif #if defined key_medusa USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA & asm_phyto2d_bal_medusa USE asmpco2bal, ONLY: & ! pCO2 balancing for MEDUSA & asm_pco2_bal 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 asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC & asm_phyto2d_bal_hadocc USE asmpco2bal, ONLY: & ! pCO2 balancing for HadOCC & asm_pco2_bal 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 #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_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_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 :: 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 # 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 ) 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 ( ( ln_phytobal ).AND. & & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. & & ( .NOT. ln_slchlnoninc ).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_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 ) ) 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. & & ( 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_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_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 ) IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .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) ) #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) ) #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) ) #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) ) #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_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 #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) ) #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) ) #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) ) #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) ) #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 IF(lwp .AND. lflush) CALL flush(numout) END SUBROUTINE asm_bgc_bal_wri !!=========================================================================== !!=========================================================================== !!=========================================================================== SUBROUTINE asm_bgc_bkg_wri( kt, knum ) !!------------------------------------------------------------------------ !! *** 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 !! !!------------------------------------------------------------------------ #if defined key_hadocc CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg' , pgrow_avg ) CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg' , ploss_avg ) CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg' , phyt_avg ) CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max' , mld_max ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut' , trn(:,:,:,jp_had_nut) ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy' , trn(:,:,:,jp_had_phy) ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo' , trn(:,:,:,jp_had_zoo) ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn' , trn(:,:,:,jp_had_pdn) ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic' , trn(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk' , trn(:,:,:,jp_had_alk) ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl' , HADOCC_CHL(:,:,:) ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:) ) #elif defined key_medusa CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg' , pgrow_avg ) CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg' , ploss_avg ) CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg' , phyt_avg ) CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max' , mld_max ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn' , trn(:,:,:,jpchn) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd' , trn(:,:,:,jpchd) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn' , trn(:,:,:,jpphn) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd' , trn(:,:,:,jpphd) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds' , trn(:,:,:,jppds) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi' , trn(:,:,:,jpzmi) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme' , trn(:,:,:,jpzme) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din' , trn(:,:,:,jpdin) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil' , trn(:,:,:,jpsil) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer' , trn(:,:,:,jpfer) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det' , trn(:,:,:,jpdet) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc' , trn(:,:,:,jpdtc) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic' , trn(:,:,:,jpdic) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk' , trn(:,:,:,jpalk) ) CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy' , trn(:,:,:,jpoxy) ) #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 #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) #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 ! Diatom chlorophyll IF ( ln_slchldiainc ) THEN zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd) 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 ! 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 ) 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(:,:) CALL ctl_stop( ' hmld_tref mixed layer requested for phyto2d assimilation,', & & ' but is not available in this version' ) 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 ) #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 #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 #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF ! IF(lwp .AND. lflush) CALL flush(numout) 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) :: zrat_phn_chn ! jpphn:jpchn ratio REAL(wp) :: zrat_phd_chd ! jpphd:jpchd ratio REAL(wp) :: zrat_pds_chd ! jppds:jpchd 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(:,:,:) #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(:,:,:) #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 #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 #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF ! IF(lwp .AND. lflush) CALL flush(numout) 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) ) #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 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(:,:) CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', & & ' but is not available in this version' ) 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 #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 #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF ! IF(lwp .AND. lflush) CALL flush(numout) 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 ! IF(lwp .AND. lflush) CALL flush(numout) 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 #if defined key_hadocc it = jp_had_nut #elif defined key_medusa it = jpdin #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 it = jpsil 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 #if defined key_hadocc it = jp_had_dic #elif defined key_medusa it = jpdic #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 #if defined key_hadocc it = jp_had_alk #elif defined key_medusa it = jpalk #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 it = jpoxy 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 #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 #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 #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 #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 #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 #if defined key_hadocc CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', & & ' Background state is taken from model rather than background file' ) #elif defined key_medusa CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', & & ' Background state is taken from model rather than background file' ) #endif 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 #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 #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 #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 #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 #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 ! IF(lwp .AND. lflush) CALL flush(numout) END SUBROUTINE bgc3d_asm_inc !!=========================================================================== END MODULE asmbgc