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_foam_medusa' : MEDUSA extras for FOAM OBS and ASM !! 'key_roam' : MEDUSA extras for carbonate cycle !! 'key_karaml' : Kara mixed layer depth !!---------------------------------------------------------------------- !! asm_bgc_check_options : check biogeochemical assimilation options !! asm_bgc_init_incs : initialise bgc increments !! asm_bgc_init_bkg : initialise bgc background !! asm_bgc_bal_wri : write out balancing increments !! asm_bgc_bkg_wri : write out bgc background !! slchltot_asm_inc : apply the slchltot increment !! spco2_asm_inc : apply the pCO2/fCO2 increment !!---------------------------------------------------------------------- USE par_kind, ONLY: wp ! kind parameters USE par_oce, ONLY: jpi, jpj, jpk ! domain array sizes USE iom ! i/o USE zdfmxl, ONLY : & ! & hmld_tref, & #if defined key_karaml & hmld_kara, & & ln_kara, & #endif & hmld, & & hmlp, & & hmlpt #if defined key_top USE trc, ONLY: & & trn, & & trb USE par_trc, ONLY: & & jptra #endif #if defined key_medusa && defined key_foam_medusa USE asmlogchlbal_medusa, ONLY: & & asm_logchl_bal_medusa USE asmpco2bal, ONLY: & & asm_pco2_bal USE par_medusa USE mocsy_mainmod ! f2pCO2 USE sms_medusa, ONLY: pgrow_avg, & & ploss_avg, & & phyt_avg, & & mld_max #elif defined key_hadocc USE asmlogchlbal_hadocc, ONLY: & & asm_logchl_bal_hadocc USE asmpco2bal, ONLY: & & asm_pco2_bal USE par_hadocc USE trc, ONLY: pgrow_avg, & & ploss_avg, & & phyt_avg, & & mld_max, & & HADOCC_CHL USE had_bgc_const, ONLY: cchl_p #endif USE asmpar ! c_asmbkg, c_asmbal, nitbkg_r, nitdin_r, nitiaustr_r, nitiaufin_r USE oce ! tsn IMPLICIT NONE PRIVATE PUBLIC asm_bgc_check_options PUBLIC asm_bgc_init_incs PUBLIC asm_bgc_init_bkg PUBLIC asm_bgc_bal_wri PUBLIC asm_bgc_bkg_wri PUBLIC slchltot_asm_inc PUBLIC spco2_asm_inc LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of balancing incs LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No slchltot increment LOGICAL, PUBLIC :: ln_slchltotbal = .FALSE. !: No slchltot balancing LOGICAL, PUBLIC :: ln_spco2inc = .FALSE. !: No pCO2 increment LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE. !: No fCO2 increment REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc !: Increment to background slchltot REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc !: Increment to background pCO2 #if defined key_top REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: slchltot_balinc !: Increment to BGC variables from logchl assim REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: spco2_balinc !: Increment to BGC variables from pCO2 assim #endif #if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg !: Background phyto growth for logchl balancing REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg !: Background phyto loss for logchl balancing REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg !: Background phyto for logchl balancing REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg !: Background max MLD for logchl balancing REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg !: Background tracer state variables #endif #if defined key_hadocc REAL(wp), DIMENSION(:,:), ALLOCATABLE :: chl_bkg !: Background surface chlorophyll REAL(wp), DIMENSION(:,:), ALLOCATABLE :: cchl_p_bkg !: Background surface carbon:chlorophyll #endif REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll increment from logchl assimilation !: <= 0 implies no maximum applied (switch turned off) !: > 0 implies maximum absolute chl increment capped at this value INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld criteria to use for biogeochemistry 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] 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 ( ( ln_balwri ).AND.( .NOT. ln_slchltotinc ).AND.( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', & & ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.') ln_balwri = .FALSE. ENDIF IF ( ( ln_slchltotbal ).AND.( .NOT. ln_slchltotinc ) ) THEN CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', & & ' Not assimilating logchl, so ln_slchltotbal will be set to .false.') ln_slchltotbal = .FALSE. ENDIF IF ( ( ln_spco2inc ).AND.( ln_sfco2inc ) ) THEN CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' ) 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 !! !!--------------------------------------------------------------------------- IF ( ln_slchltotinc ) THEN ALLOCATE( slchltot_bkginc(jpi,jpj) ) slchltot_bkginc(:,:) = 0.0 #if defined key_top ALLOCATE( slchltot_balinc(jpi,jpj,jpk,jptra) ) slchltot_balinc(:,:,:,:) = 0.0 #endif ENDIF IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN ALLOCATE( spco2_bkginc(jpi,jpj) ) spco2_bkginc(:,:) = 0.0 #if defined key_top ALLOCATE( spco2_balinc(jpi,jpj,jpk,jptra) ) spco2_balinc(:,:,:,:) = 0.0 #endif ENDIF IF ( ln_slchltotinc ) THEN CALL iom_get( knum, jpdom_autoglo, 'bckinslchltot', slchltot_bkginc(:,:), 1 ) ! Apply the masks slchltot_bkginc(:,:) = slchltot_bkginc(:,:) * tmask(:,:,1) ! Set missing increments to 0.0 rather than 1e+20 ! to allow for differences in masks WHERE( ABS( slchltot_bkginc(:,:) ) > 1.0e+10 ) slchltot_bkginc(:,:) = 0.0 ENDIF IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN IF ( ln_spco2inc ) THEN CALL iom_get( knum, jpdom_autoglo, 'bckinspco2', spco2_bkginc(:,:), 1 ) ELSE IF ( ln_sfco2inc ) THEN CALL iom_get( knum, jpdom_autoglo, 'bckinsfco2', spco2_bkginc(:,:), 1 ) ENDIF ! Apply the masks spco2_bkginc(:,:) = spco2_bkginc(:,:) * tmask(:,:,1) ! Set missing increments to 0.0 rather than 1e+20 ! to allow for differences in masks WHERE( ABS( spco2_bkginc(:,:) ) > 1.0e+10 ) spco2_bkginc(:,:) = 0.0 ENDIF END SUBROUTINE asm_bgc_init_incs !!====================================================================== !!====================================================================== !!====================================================================== 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_hadocc || (defined key_medusa && defined key_foam_medusa) IF ( ln_slchltotinc ) 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) ) ALLOCATE( cchl_p_bkg(jpi,jpj) ) chl_bkg(:,:) = 0.0 cchl_p_bkg(:,:) = 0.0 #endif !-------------------------------------------------------------------- ! Read background variables for slchltot 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(:,:,1) cchl_p_bkg(:,:) = cchl_p_bkg(:,:) * tmask(:,:,1) #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) ) #endif IF ( ln_slchltotbal ) 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_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) ) 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 ) 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 ) 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 #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 ) THEN #if defined key_medusa CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', slchltot_balinc(:,:,:,jpchn) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', slchltot_balinc(:,:,:,jpchd) ) IF ( ln_slchltotbal ) THEN CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', slchltot_balinc(:,:,:,jpphn) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', slchltot_balinc(:,:,:,jpphd) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', slchltot_balinc(:,:,:,jppds) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', slchltot_balinc(:,:,:,jpzmi) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', slchltot_balinc(:,:,:,jpzme) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', slchltot_balinc(:,:,:,jpdin) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', slchltot_balinc(:,:,:,jpsil) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', slchltot_balinc(:,:,:,jpfer) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', slchltot_balinc(:,:,:,jpdet) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', slchltot_balinc(:,:,:,jpdtc) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', slchltot_balinc(:,:,:,jpdic) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', slchltot_balinc(:,:,:,jpalk) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', slchltot_balinc(:,:,:,jpoxy) ) ENDIF #elif defined key_hadocc CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', slchltot_balinc(:,:,:,jp_had_phy) ) IF ( ln_slchltotbal ) THEN CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', slchltot_balinc(:,:,:,jp_had_nut) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', slchltot_balinc(:,:,:,jp_had_zoo) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', slchltot_balinc(:,:,:,jp_had_pdn) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', slchltot_balinc(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', slchltot_balinc(:,:,:,jp_had_alk) ) ENDIF #endif ENDIF IF ( ln_spco2inc ) THEN #if defined key_medusa CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', spco2_balinc(:,:,:,jpdic) ) CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', spco2_balinc(:,:,:,jpalk) ) #elif defined key_hadocc CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', spco2_balinc(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', spco2_balinc(:,:,:,jp_had_alk) ) #endif ELSE IF ( ln_sfco2inc ) THEN #if defined key_medusa CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', spco2_balinc(:,:,:,jpdic) ) CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', spco2_balinc(:,:,:,jpalk) ) #elif defined key_hadocc CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', spco2_balinc(:,:,:,jp_had_dic) ) CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', spco2_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 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(:,:,1) ) CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,1) ) #elif defined key_medusa && defined key_foam_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 slchltot_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau ) !!---------------------------------------------------------------------- !! *** ROUTINE logchl_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 :: jk ! Loop counter INTEGER :: it ! Index REAL(wp) :: zincwgt ! IAU weight for current time step REAL(wp) :: zincper ! IAU interval in seconds !!---------------------------------------------------------------------- IF ( kt <= nit000 ) THEN zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt #if defined key_medusa && defined key_foam_medusa CALL asm_logchl_bal_medusa( slchltot_bkginc, zincper, mld_choice_bgc, & & rn_maxchlinc, ln_slchltotbal, ll_asmdin, & & pgrow_avg_bkg, ploss_avg_bkg, & & phyt_avg_bkg, mld_max_bkg, & & tracer_bkg, slchltot_balinc ) #elif defined key_hadocc CALL asm_logchl_bal_hadocc( slchltot_bkginc, zincper, mld_choice_bgc, & & rn_maxchlinc, ln_slchltotbal, ll_asmdin, & & pgrow_avg_bkg, ploss_avg_bkg, & & phyt_avg_bkg, mld_max_bkg, & & chl_bkg, cchl_p_bkg, & & tracer_bkg, slchltot_balinc ) #else CALL ctl_stop( 'Attempting to assimilate slchltot, ', & & '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,*) 'logchl_asm_inc : logchl 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 && defined key_foam_medusa DO jk = 1, jpkm1 trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & & slchltot_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & & slchltot_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt END DO #elif defined key_hadocc DO jk = 1, jpkm1 trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & & slchltot_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & & slchltot_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt END DO #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 #if defined key_medusa && defined key_foam_medusa ! 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 chlorophyll assimilation', & & ' Background state is taken from model rather than background file' ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & slchltot_balinc(:,:,:,jp_msa0:jp_msa1) trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) #elif defined key_hadocc ! Initialize the now fields with the background + increment ! Background currently is what the model is initialised with CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', & & ' Background state is taken from model rather than background file' ) trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & & slchltot_balinc(:,:,:,jp_had0:jp_had1) trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF ! END SUBROUTINE slchltot_asm_inc !!====================================================================== !!====================================================================== !!====================================================================== SUBROUTINE spco2_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( spco2_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(:,:) = spco2_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 logchl balancing if required IF ( ln_slchltotinc .AND. ln_slchltotbal ) THEN dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + slchltot_balinc(:,:,1,jpdic) alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + slchltot_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(:,:), & & spco2_balinc(:,:,1,jpdic), spco2_balinc(:,:,1,jpalk) ) #elif defined key_hadocc ! Account for slchltot balancing if required IF ( ln_slchltotinc .AND. ln_slchltotbal ) THEN dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + slchltot_balinc(:,:,1,jp_had_dic) alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + slchltot_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(:,:), & & spco2_balinc(:,:,1,jp_had_dic), spco2_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 && defined key_foam_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 spco2_balinc(ji,jj,jk,:) = spco2_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 && defined key_foam_medusa DO jk = 1, jpkm1 trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & & spco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & & spco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt END DO #elif defined key_hadocc DO jk = 1, jpkm1 trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & & spco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & & spco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt END DO #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 #if defined key_medusa && defined key_foam_medusa ! 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 pCO2 assimilation', & & ' Background state is taken from model rather than background file' ) trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & & spco2_balinc(:,:,:,jp_msa0:jp_msa1) trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) #elif defined key_hadocc ! Initialize the now fields with the background + increment ! Background currently is what the model is initialised with CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', & & ' Background state is taken from model rather than background file' ) trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & & spco2_balinc(:,:,:,jp_had0:jp_had1) trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) #endif ! Do not deallocate arrays - needed by asm_bgc_bal_wri ! which is called at end of model run ENDIF ! ENDIF ! END SUBROUTINE spco2_asm_inc !!====================================================================== END MODULE asmbgc