Changeset 9296


Ignore:
Timestamp:
2018-01-31T16:27:16+01:00 (2 years ago)
Author:
dford
Message:

Move most of the BGC assimilation code into its own module, to make it more portable across versions/configurations. Partially apply the new naming conventions for BGC obs.

Location:
branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM
Files:
1 added
1 deleted
5 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref

    r9292 r9296  
    12941294    ln_asmiau = .false.    !  Logical switch for Incremental Analysis Updating (IAU) 
    12951295    ln_seaiceinc = .false. !  Logical switch for applying sea ice increments 
    1296     ln_logchlinc = .false. !  Logical switch for applying logchl increments 
    1297     ln_logchlbal = .false. !  Logical switch for applying logchl multivariate balancing 
    1298     ln_fco2inc = .false.   !  Logical switch for applying fCO2 increments 
    1299     ln_pco2inc = .false.   !  Logical switch for applying pCO2 increments 
     1296    ln_slchltotinc = .false. !  Logical switch for applying slchltot increments 
     1297    ln_slchltotbal = .false. !  Logical switch for applying slchltot multivariate balancing 
     1298    ln_sfco2inc = .false.   !  Logical switch for applying sfCO2 increments 
     1299    ln_spco2inc = .false.   !  Logical switch for applying spCO2 increments 
    13001300    ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing 
    13011301    nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1] 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r9292 r9296  
    5050   USE ice 
    5151#endif 
    52 #if defined key_hadocc 
    53    USE trc, ONLY: trn, & 
    54       &     pgrow_avg, & 
    55       &     ploss_avg, & 
    56       &     phyt_avg,  & 
    57       &     mld_max,   & 
    58       &     HADOCC_CHL 
    59    USE had_bgc_const, ONLY: cchl_p 
    60    USE par_hadocc 
    61 #elif defined key_medusa && defined key_foam_medusa 
    62    USE trc, ONLY: trn 
    63    USE sms_medusa, ONLY: pgrow_avg, & 
    64       &                  ploss_avg, & 
    65       &                  phyt_avg,  & 
    66       &                  mld_max 
    67    USE par_medusa 
     52#if defined key_top 
     53   USE asmbgc, ONLY: asm_bgc_bkg_wri 
    6854#endif 
    6955   IMPLICIT NONE 
     
    138124!            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
    139125            CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               ) 
    140 #if defined key_hadocc 
    141             CALL iom_rstput( kt, nitbkg_r, inum, 'pgrow_avg'   , pgrow_avg             ) 
    142             CALL iom_rstput( kt, nitbkg_r, inum, 'ploss_avg'   , ploss_avg             ) 
    143             CALL iom_rstput( kt, nitbkg_r, inum, 'phyt_avg'    , phyt_avg              ) 
    144             CALL iom_rstput( kt, nitbkg_r, inum, 'mld_max'     , mld_max               ) 
    145             CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) ) 
    146             CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) ) 
    147             CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_zoo'  , trn(:,:,:,jp_had_zoo) ) 
    148             CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_pdn'  , trn(:,:,:,jp_had_pdn) ) 
    149             CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_dic'  , trn(:,:,:,jp_had_dic) ) 
    150             CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_alk'  , trn(:,:,:,jp_had_alk) ) 
    151             CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_chl'  , HADOCC_CHL(:,:,1)     ) 
    152             CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_cchl' , cchl_p(:,:,1)         ) 
    153 #elif defined key_medusa && defined key_foam_medusa 
    154             CALL iom_rstput( kt, nitbkg_r, inum, 'pgrow_avg'   , pgrow_avg        ) 
    155             CALL iom_rstput( kt, nitbkg_r, inum, 'ploss_avg'   , ploss_avg        ) 
    156             CALL iom_rstput( kt, nitbkg_r, inum, 'phyt_avg'    , phyt_avg         ) 
    157             CALL iom_rstput( kt, nitbkg_r, inum, 'mld_max'     , mld_max          ) 
    158             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_chn'  , trn(:,:,:,jpchn) ) 
    159             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_chd'  , trn(:,:,:,jpchd) ) 
    160             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_phn'  , trn(:,:,:,jpphn) ) 
    161             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_phd'  , trn(:,:,:,jpphd) ) 
    162             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_pds'  , trn(:,:,:,jppds) ) 
    163             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_zmi'  , trn(:,:,:,jpzmi) ) 
    164             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_zme'  , trn(:,:,:,jpzme) ) 
    165             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_din'  , trn(:,:,:,jpdin) ) 
    166             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_sil'  , trn(:,:,:,jpsil) ) 
    167             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_fer'  , trn(:,:,:,jpfer) ) 
    168             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_det'  , trn(:,:,:,jpdet) ) 
    169             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_dtc'  , trn(:,:,:,jpdtc) ) 
    170             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_dic'  , trn(:,:,:,jpdic) ) 
    171             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_alk'  , trn(:,:,:,jpalk) ) 
    172             CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_oxy'  , trn(:,:,:,jpoxy) ) 
     126            ! 
     127#if defined key_top 
     128            CALL asm_bgc_bkg_wri( kt, inum ) 
     129            ! 
    173130#endif 
    174             ! 
    175131            CALL iom_close( inum ) 
    176132         ENDIF 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r9292 r9296  
    2222   !!   ssh_asm_inc    : Apply the SSH increment 
    2323   !!   seaice_asm_inc : Apply the seaice increment 
    24    !!   logchl_asm_inc : Apply the logchl increment 
    25    !!   pco2_asm_inc   : Apply the pCO2/fCO2 increment 
    2624   !!---------------------------------------------------------------------- 
    2725   USE wrk_nemo         ! Memory Allocation 
     
    4644#endif 
    4745   USE sbc_oce          ! Surface boundary condition variables. 
    48    USE zdfmxl, ONLY :  &   
    49 !   &  hmld_tref,       &    
    50 #if defined key_karaml 
    51    &  hmld_kara,       & 
    52    &  ln_kara,         & 
    53 #endif    
    54    &  hmld,            &  
    55    &  hmlp,            & 
    56    &  hmlpt  
    57 #if defined key_top 
    58    USE trc, ONLY: & 
    59       & trn,      & 
    60       & trb 
    61    USE par_trc, ONLY: & 
    62       & jptra 
    63 #endif 
    64 #if defined key_medusa && defined key_foam_medusa 
    65    USE asmlogchlbal_medusa, ONLY: & 
    66       & asm_logchl_bal_medusa 
    67    USE asmpco2bal, ONLY: & 
    68       & asm_pco2_bal 
    69    USE par_medusa 
    70    USE mocsy_mainmod 
    71 #elif defined key_hadocc 
    72    USE asmlogchlbal_hadocc, ONLY: & 
    73       & asm_logchl_bal_hadocc 
    74    USE asmpco2bal, ONLY: & 
    75       & asm_pco2_bal 
    76    USE par_hadocc 
    77 #endif 
     46   USE asmbgc           ! Biogeochemistry assimilation 
    7847 
    7948   IMPLICIT NONE 
     
    8655   PUBLIC   ssh_asm_inc    !: Apply the SSH increment 
    8756   PUBLIC   seaice_asm_inc !: Apply the seaice increment 
    88    PUBLIC   logchl_asm_inc !: Apply the logchl increment 
    89    PUBLIC   pco2_asm_inc   !: Apply the pCO2/fCO2 increment 
     57   PUBLIC   bgc_asm_inc    !: Apply the biogeochemistry increments 
    9058 
    9159#if defined key_asminc 
     
    9563#endif 
    9664   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields 
    97    LOGICAL, PUBLIC :: ln_balwri = .FALSE.      !: No output of the assimilation balancing increments 
    9865   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment 
    9966   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization 
     
    10269   LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment 
    10370   LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment 
    104    LOGICAL, PUBLIC :: ln_logchlinc = .FALSE.   !: No log10(chlorophyll) increment 
    105    LOGICAL, PUBLIC :: ln_logchlbal = .FALSE.   !: Don't apply multivariate log10(chlorophyll) balancing 
    106    LOGICAL, PUBLIC :: ln_pco2inc = .FALSE.     !: No pCO2 increment 
    107    LOGICAL, PUBLIC :: ln_fco2inc = .FALSE.     !: No fCO2 increment 
     71   LOGICAL, PUBLIC :: lk_bgcinc = .FALSE.      !: No biogeochemistry increments 
    10872   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check 
    10973   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
     
    13094   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    13195   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
    132  
    133    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl 
    134    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pco2_bkginc   !: Increment to background pCO2 
    135 #if defined key_top 
    136    REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim 
    137    REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc    !: Increment to BGC variables from pCO2 assim 
    138 #endif 
    139 #if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 
    140    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg  !: Background phyto growth for logchl balancing 
    141    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg  !: Background phyto loss for logchl balancing 
    142    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg   !: Background phyto for logchl balancing 
    143    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg    !: Background max MLD for logchl balancing 
    144    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg     !: Background tracer state variables 
    145 #endif 
    146 #if defined key_hadocc 
    147    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: chl_bkg        !: Background surface chlorophyll 
    148    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: cchl_p_bkg     !: Background surface carbon:chlorophyll 
    149 #endif 
    150    REAL(wp) :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation 
    151                                       !: <= 0 implies no maximum applied (switch turned off) 
    152                                       !:  > 0 implies maximum absolute chl increment capped at this value 
    153  
    154    INTEGER :: mld_choice_bgc    = 1   !: choice of mld criteria to use for biogeochemistry assimilation 
    155                                       !: 1) hmld      - Turbocline/mixing depth                           [W points] 
    156                                       !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points] 
    157                                       !: 3) hmld_kara - Kara MLD                                          [Interpolated] 
    158                                       !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 
    159                                       !: 5) hmlpt     - Density criterion (0.01 kg/m^3 change from 10m)   [T points] 
    16096 
    16197   !! * Substitutions 
     
    202138      NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri,                           & 
    203139         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    204          &                 ln_logchlinc, ln_logchlbal,                     & 
    205          &                 ln_pco2inc, ln_fco2inc,                         & 
     140         &                 ln_slchltotinc, ln_slchltotbal,                 & 
     141         &                 ln_spco2inc, ln_sfco2inc,                       & 
    206142         &                 ln_asmdin, ln_asmiau,                           & 
    207143         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
     
    239175         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin 
    240176         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc 
    241          WRITE(numout,*) '      Logical switch for applying logchl increments         ln_logchlinc = ', ln_logchlinc 
    242          WRITE(numout,*) '      Logical switch for logchl multivariate balancing      ln_logchlbal = ', ln_logchlbal 
    243          WRITE(numout,*) '      Logical switch for applying pCO2 increments             ln_pco2inc = ', ln_pco2inc 
    244          WRITE(numout,*) '      Logical switch for applying fCO2 increments             ln_fco2inc = ', ln_fco2inc 
     177         WRITE(numout,*) '      Logical switch for applying slchltot increments     ln_slchltotinc = ', ln_slchltotinc 
     178         WRITE(numout,*) '      Logical switch for slchltot multivariate balancing  ln_slchltotbal = ', ln_slchltotbal 
     179         WRITE(numout,*) '      Logical switch for applying pCO2 increments            ln_spco2inc = ', ln_spco2inc 
     180         WRITE(numout,*) '      Logical switch for applying fCO2 increments            ln_sfco2inc = ', ln_sfco2inc 
    245181         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau 
    246182         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg 
     
    289225         WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date 
    290226      ENDIF 
     227       
     228      IF ( ln_slchltotinc .OR. ln_spco2inc .OR. ln_sfco2inc ) THEN 
     229         lk_bgcinc = .TRUE. 
     230      ENDIF 
    291231 
    292232      IF ( nacc /= 0 ) & 
     
    301241      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 
    302242         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. & 
    303          &        ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) )) & 
     243         &        ( lk_bgcinc ) )) & 
    304244         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
    305          &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc is set to .true.', & 
     245         &                ' ln_(bgc-variable)inc is set to .true.', & 
    306246         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', & 
    307247         &                ' Inconsistent options') 
     
    312252 
    313253      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 
    314          & .AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) )  & 
     254         & .AND.( .NOT. lk_bgcinc ) )  & 
    315255         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
    316          &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc are set to .false. :', & 
     256         &                ' ln_(bgc-variable)inc are set to .false. :', & 
    317257         &                ' The assimilation increments are not applied') 
    318258 
     
    334274         &                ' Background time step for Direct Initialization is outside', & 
    335275         &                ' the cycle interval') 
    336  
    337       IF ( ( ln_balwri ).AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) ) THEN 
    338          CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', & 
    339             &           ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.') 
    340          ln_balwri = .FALSE. 
    341       ENDIF 
    342  
    343       IF ( ( ln_logchlbal ).AND.( .NOT. ln_logchlinc ) ) THEN 
    344          CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', & 
    345             &           ' Not assimilating logchl, so ln_logchlbal will be set to .false.') 
    346          ln_logchlbal = .FALSE. 
    347       ENDIF 
    348  
    349       IF ( ( ln_pco2inc ).AND.( ln_fco2inc ) ) THEN 
    350          CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' ) 
    351       ENDIF 
     276       
     277      IF ( lk_bgcinc ) CALL asm_bgc_check_options 
    352278 
    353279      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further 
     
    446372      ssh_iau(:,:)    = 0.0 
    447373#endif 
    448       IF ( ln_logchlinc ) THEN 
    449          ALLOCATE( logchl_bkginc(jpi,jpj) ) 
    450          logchl_bkginc(:,:) = 0.0 
    451 #if defined key_top 
    452          ALLOCATE( logchl_balinc(jpi,jpj,jpk,jptra) ) 
    453          logchl_balinc(:,:,:,:) = 0.0 
    454 #endif 
    455       ENDIF 
    456       IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 
    457          ALLOCATE( pco2_bkginc(jpi,jpj) ) 
    458          pco2_bkginc(:,:) = 0.0 
    459 #if defined key_top 
    460          ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) ) 
    461          pco2_balinc(:,:,:,:) = 0.0 
    462 #endif 
    463       ENDIF 
    464374      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 
    465          &  .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN 
     375         &  .OR.( lk_bgcinc ) ) THEN 
    466376 
    467377         !-------------------------------------------------------------------- 
     
    537447         ENDIF 
    538448 
    539          IF ( ln_logchlinc ) THEN 
    540             CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl', logchl_bkginc(:,:), 1 ) 
    541             ! Apply the masks 
    542             logchl_bkginc(:,:) = logchl_bkginc(:,:) * tmask(:,:,1) 
    543             ! Set missing increments to 0.0 rather than 1e+20 
    544             ! to allow for differences in masks 
    545             WHERE( ABS( logchl_bkginc(:,:) ) > 1.0e+10 ) logchl_bkginc(:,:) = 0.0 
    546          ENDIF 
    547  
    548          IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 
    549             IF ( ln_pco2inc ) THEN 
    550                CALL iom_get( inum, jpdom_autoglo, 'bckinpco2', pco2_bkginc(:,:), 1 ) 
    551             ELSE IF ( ln_fco2inc ) THEN 
    552                CALL iom_get( inum, jpdom_autoglo, 'bckinfco2', pco2_bkginc(:,:), 1 ) 
    553             ENDIF 
    554             ! Apply the masks 
    555             pco2_bkginc(:,:) = pco2_bkginc(:,:) * tmask(:,:,1) 
    556             ! Set missing increments to 0.0 rather than 1e+20 
    557             ! to allow for differences in masks 
    558             WHERE( ABS( pco2_bkginc(:,:) ) > 1.0e+10 ) pco2_bkginc(:,:) = 0.0 
     449         IF ( lk_bgcinc ) THEN 
     450            CALL asm_bgc_init_incs( inum ) 
    559451         ENDIF 
    560452 
     
    670562 
    671563      ENDIF 
    672  
    673 #if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 
    674       IF ( ln_logchlinc ) THEN 
    675  
    676          ALLOCATE( pgrow_avg_bkg(jpi,jpj)        ) 
    677          ALLOCATE( ploss_avg_bkg(jpi,jpj)        ) 
    678          ALLOCATE( phyt_avg_bkg(jpi,jpj)         ) 
    679          ALLOCATE( mld_max_bkg(jpi,jpj)          ) 
    680          ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) 
    681          pgrow_avg_bkg(:,:)  = 0.0 
    682          ploss_avg_bkg(:,:)  = 0.0 
    683          phyt_avg_bkg(:,:)   = 0.0 
    684          mld_max_bkg(:,:)    = 0.0 
    685          tracer_bkg(:,:,:,:) = 0.0 
    686  
    687 #if defined key_hadocc 
    688          ALLOCATE( chl_bkg(jpi,jpj)    ) 
    689          ALLOCATE( cchl_p_bkg(jpi,jpj) ) 
    690          chl_bkg(:,:)    = 0.0 
    691          cchl_p_bkg(:,:) = 0.0 
    692 #endif 
    693564          
    694          !-------------------------------------------------------------------- 
    695          ! Read background variables for logchl assimilation 
    696          ! Some only required if performing balancing 
    697          !-------------------------------------------------------------------- 
    698  
    699          CALL iom_open( c_asmbkg, inum ) 
    700  
    701 #if defined key_hadocc 
    702          CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl',  chl_bkg    ) 
    703          CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg ) 
    704          chl_bkg(:,:)    = chl_bkg(:,:)    * tmask(:,:,1) 
    705          cchl_p_bkg(:,:) = cchl_p_bkg(:,:) * tmask(:,:,1) 
    706 #elif defined key_medusa 
    707          CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) ) 
    708          CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) ) 
    709 #endif 
    710           
    711          IF ( ln_logchlbal ) THEN 
    712  
    713             CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg ) 
    714             CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg ) 
    715             CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  ) 
    716             CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   ) 
    717             pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1) 
    718             ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1) 
    719             phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1) 
    720             mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1) 
    721  
    722 #if defined key_hadocc 
    723             CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) ) 
    724             CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) ) 
    725             CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) ) 
    726             CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) ) 
    727             CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 
    728             CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) 
    729 #elif defined key_medusa 
    730             CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) ) 
    731             CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) ) 
    732             CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) ) 
    733             CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) ) 
    734             CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) ) 
    735             CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) ) 
    736             CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) ) 
    737             CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) ) 
    738             CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) ) 
    739             CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) ) 
    740             CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 
    741             CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 
    742             CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) ) 
    743 #endif 
    744          ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 
    745 #if defined key_hadocc 
    746             CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 
    747             CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) 
    748 #elif defined key_medusa 
    749             CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 
    750             CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 
    751 #endif 
    752             CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   ) 
    753             mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) 
    754          ENDIF 
    755  
    756          CALL iom_close( inum ) 
    757           
    758          DO jt = 1, jptra 
    759             tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 
    760          END DO 
    761        
    762       ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 
    763  
    764          ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) 
    765          ALLOCATE( mld_max_bkg(jpi,jpj)          ) 
    766          tracer_bkg(:,:,:,:) = 0.0 
    767          mld_max_bkg(:,:)    = 0.0 
    768  
    769          CALL iom_open( c_asmbkg, inum ) 
    770           
    771 #if defined key_hadocc 
    772          CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 
    773          CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) 
    774 #elif defined key_medusa 
    775          CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 
    776          CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 
    777 #endif 
    778          CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   ) 
    779  
    780          CALL iom_close( inum ) 
    781           
    782          DO jt = 1, jptra 
    783             tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 
    784          END DO 
    785          mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) 
    786        
    787       ENDIF 
    788 #endif 
     565      IF ( lk_bgcinc ) THEN 
     566         CALL asm_bgc_init_bkg 
     567      ENDIF 
    789568      ! 
    790569   END SUBROUTINE asm_inc_init 
     
    14081187   END SUBROUTINE seaice_asm_inc 
    14091188 
    1410    SUBROUTINE logchl_asm_inc( kt ) 
    1411       !!---------------------------------------------------------------------- 
    1412       !!                    ***  ROUTINE logchl_asm_inc  *** 
     1189 
     1190   SUBROUTINE bgc_asm_inc( kt ) 
     1191      !!---------------------------------------------------------------------- 
     1192      !!                    ***  ROUTINE bgc_asm_inc  *** 
    14131193      !!           
    1414       !! ** Purpose : Apply the chlorophyll assimilation increments. 
    1415       !! 
    1416       !! ** Method  : Calculate increments to state variables using nitrogen 
    1417       !!              balancing. 
    1418       !!              Direct initialization or Incremental Analysis Updating. 
    1419       !! 
    1420       !! ** Action  :  
    1421       !!---------------------------------------------------------------------- 
    1422       INTEGER, INTENT(IN) :: kt   ! Current time step 
    1423       ! 
    1424       INTEGER  :: jk              ! Loop counter 
    1425       INTEGER  :: it              ! Index 
    1426       REAL(wp) :: zincwgt         ! IAU weight for current time step 
    1427       REAL(wp) :: zincper         ! IAU interval in seconds 
    1428       !!---------------------------------------------------------------------- 
    1429  
    1430       IF ( kt <= nit000 ) THEN 
    1431  
    1432          zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt 
    1433  
    1434 #if defined key_medusa && defined key_foam_medusa 
    1435          CALL asm_logchl_bal_medusa( logchl_bkginc, zincper, mld_choice_bgc, & 
    1436             &                        rn_maxchlinc, ln_logchlbal, ln_asmdin,  & 
    1437             &                        pgrow_avg_bkg, ploss_avg_bkg,           & 
    1438             &                        phyt_avg_bkg, mld_max_bkg,              & 
    1439             &                        tracer_bkg, logchl_balinc ) 
    1440 #elif defined key_hadocc 
    1441          CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, & 
    1442             &                        rn_maxchlinc, ln_logchlbal, ln_asmdin,  & 
    1443             &                        pgrow_avg_bkg, ploss_avg_bkg,           & 
    1444             &                        phyt_avg_bkg, mld_max_bkg,              & 
    1445             &                        chl_bkg, cchl_p_bkg,                    & 
    1446             &                        tracer_bkg, logchl_balinc ) 
    1447 #else 
    1448          CALL ctl_stop( 'Attempting to assimilate logchl, ', & 
    1449             &           'but not defined a biogeochemical model' ) 
    1450 #endif 
    1451  
    1452       ENDIF 
    1453  
    1454       IF ( ln_asmiau ) THEN 
    1455  
    1456          !-------------------------------------------------------------------- 
    1457          ! Incremental Analysis Updating 
    1458          !-------------------------------------------------------------------- 
    1459  
    1460          IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
    1461  
    1462             it = kt - nit000 + 1 
    1463             zincwgt = wgtiau(it)      ! IAU weight for the current time step 
    1464             ! note this is not a tendency so should not be divided by rdt 
    1465  
    1466             IF(lwp) THEN 
    1467                WRITE(numout,*)  
    1468                WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', & 
    1469                   &  kt,' with IAU weight = ', wgtiau(it) 
    1470                WRITE(numout,*) '~~~~~~~~~~~~' 
    1471             ENDIF 
    1472  
    1473             ! Update the biogeochemical variables 
    1474             ! Add directly to trn and trb, rather than to tra, because tra gets 
    1475             ! reset to zero at the start of trc_stp, called after this routine 
    1476 #if defined key_medusa && defined key_foam_medusa 
    1477             DO jk = 1, jpkm1 
    1478                trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 
    1479                   &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
    1480                trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 
    1481                   &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
    1482             END DO 
    1483 #elif defined key_hadocc 
    1484             DO jk = 1, jpkm1 
    1485                trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 
    1486                   &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
    1487                trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 
    1488                   &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
    1489             END DO 
    1490 #endif 
    1491  
    1492             ! Do not deallocate arrays - needed by asm_bal_wri 
    1493             ! which is called at end of model run 
    1494          ENDIF 
    1495  
    1496       ELSEIF ( ln_asmdin ) THEN  
    1497  
    1498          !-------------------------------------------------------------------- 
    1499          ! Direct Initialization 
    1500          !-------------------------------------------------------------------- 
    1501           
    1502          IF ( kt == nitdin_r ) THEN 
    1503  
    1504             neuler = 0                    ! Force Euler forward step 
    1505  
    1506 #if defined key_medusa && defined key_foam_medusa 
    1507             ! Initialize the now fields with the background + increment 
    1508             ! Background currently is what the model is initialised with 
    1509             CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', & 
    1510                &           ' Background state is taken from model rather than background file' ) 
    1511             trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 
    1512                &                         logchl_balinc(:,:,:,jp_msa0:jp_msa1) 
    1513             trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 
    1514 #elif defined key_hadocc 
    1515             ! Initialize the now fields with the background + increment 
    1516             ! Background currently is what the model is initialised with 
    1517             CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', & 
    1518                &           ' Background state is taken from model rather than background file' ) 
    1519             trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 
    1520                &                         logchl_balinc(:,:,:,jp_had0:jp_had1) 
    1521             trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 
    1522 #endif 
    1523   
    1524             ! Do not deallocate arrays - needed by asm_bal_wri 
    1525             ! which is called at end of model run 
    1526          ENDIF 
    1527          ! 
    1528       ENDIF 
    1529       ! 
    1530    END SUBROUTINE logchl_asm_inc 
    1531  
    1532  
    1533    SUBROUTINE pco2_asm_inc( kt ) 
    1534       !!---------------------------------------------------------------------- 
    1535       !!                    ***  ROUTINE pco2_asm_inc  *** 
    1536       !!           
    1537       !! ** Purpose : Apply the pco2/fco2 assimilation increments. 
    1538       !! 
    1539       !! ** Method  : Calculate increments to state variables using carbon 
    1540       !!              balancing. 
    1541       !!              Direct initialization or Incremental Analysis Updating. 
    1542       !! 
    1543       !! ** Action  :  
    1544       !!---------------------------------------------------------------------- 
    1545       INTEGER, INTENT(IN)                   :: kt           ! Current time step 
    1546       ! 
    1547       INTEGER                               :: ji, jj, jk   ! Loop counters 
    1548       INTEGER                               :: it           ! Index 
    1549       INTEGER                               :: jkmax        ! Index of mixed layer depth 
    1550       REAL(wp)                              :: zincwgt      ! IAU weight for current time step 
    1551       REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing 
    1552       REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion 
    1553       REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background 
    1554       REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background 
    1555       REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background 
    1556       REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background 
    1557       REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure 
    1558       REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure 
    1559  
    1560       ! Coefficients for fCO2 to pCO2 conversion 
    1561       REAL(wp)                              :: zcoef_fco2_1 = -1636.75 
    1562       REAL(wp)                              :: zcoef_fco2_2 = 12.0408 
    1563       REAL(wp)                              :: zcoef_fco2_3 = 0.0327957 
    1564       REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528 
    1565       REAL(wp)                              :: zcoef_fco2_5 = 2.0 
    1566       REAL(wp)                              :: zcoef_fco2_6 = 57.7 
    1567       REAL(wp)                              :: zcoef_fco2_7 = 0.118 
    1568       REAL(wp)                              :: zcoef_fco2_8 = 82.0578 
    1569       !!---------------------------------------------------------------------- 
    1570  
    1571       IF ( kt <= nit000 ) THEN 
    1572  
    1573          !-------------------------------------------------------------------- 
    1574          ! DIC and alkalinity balancing 
    1575          !-------------------------------------------------------------------- 
    1576  
    1577          IF ( ln_fco2inc ) THEN 
    1578 #if defined key_medusa && defined key_roam 
    1579             ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine 
    1580             patm(1) = 1.0 
    1581             phyd(1) = 0.0 
    1582             DO jj = 1, jpj 
    1583                DO ji = 1, jpi 
    1584                   CALL f2pCO2( pco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) ) 
    1585                END DO 
    1586             END DO 
    1587 #else 
    1588             ! If assimilating fCO2, then convert to pCO2 using temperature 
    1589             ! See flux_gas.F90 within HadOCC for details of calculation 
    1590             pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) /                                                             & 
    1591                &                    EXP((zcoef_fco2_1                                                            + & 
    1592                &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - & 
    1593                &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + & 
    1594                &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + & 
    1595                &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / & 
    1596                &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0))) 
    1597 #endif 
    1598          ELSE 
    1599             pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) 
    1600          ENDIF 
    1601  
    1602          ! Account for physics assimilation if required 
    1603          IF ( ln_trainc ) THEN 
    1604             tem_bkg_temp(:,:) = tsn(:,:,1,1) + t_bkginc(:,:,1) 
    1605             sal_bkg_temp(:,:) = tsn(:,:,1,2) + s_bkginc(:,:,1) 
    1606          ELSE 
    1607             tem_bkg_temp(:,:) = tsn(:,:,1,1) 
    1608             sal_bkg_temp(:,:) = tsn(:,:,1,2) 
    1609          ENDIF 
    1610  
    1611 #if defined key_medusa 
    1612          ! Account for logchl balancing if required 
    1613          IF ( ln_logchlinc .AND. ln_logchlbal ) THEN 
    1614             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + logchl_balinc(:,:,1,jpdic) 
    1615             alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + logchl_balinc(:,:,1,jpalk) 
    1616          ELSE 
    1617             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) 
    1618             alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) 
    1619          ENDIF 
    1620  
    1621          CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 
    1622             &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        & 
    1623             &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) ) 
    1624  
    1625 #elif defined key_hadocc 
    1626          ! Account for logchl balancing if required 
    1627          IF ( ln_logchlinc .AND. ln_logchlbal ) THEN 
    1628             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + logchl_balinc(:,:,1,jp_had_dic) 
    1629             alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + logchl_balinc(:,:,1,jp_had_alk) 
    1630          ELSE 
    1631             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) 
    1632             alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) 
    1633          ENDIF 
    1634  
    1635          CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 
    1636             &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        & 
    1637             &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) ) 
    1638  
    1639 #else 
    1640          CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', & 
    1641             &           'but not defined a biogeochemical model' ) 
    1642 #endif 
    1643  
    1644          ! Select mixed layer 
    1645          IF ( ln_asmdin ) THEN 
    1646 #if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 
    1647             CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', & 
    1648                &           ' Mixed layer depth taken to be background maximum mld_max_bkg' ) 
    1649             zmld(:,:) = mld_max_bkg(:,:) 
    1650 #else 
    1651             CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' ) 
    1652 #endif 
    1653          ELSE 
    1654             SELECT CASE( mld_choice_bgc ) 
    1655             CASE ( 1 )                   ! Turbocline/mixing depth [W points] 
    1656                zmld(:,:) = hmld(:,:) 
    1657             CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points] 
    1658                zmld(:,:) = hmlp(:,:) 
    1659             CASE ( 3 )                   ! Kara MLD [Interpolated] 
    1660 #if defined key_karaml 
    1661                IF ( ln_kara ) THEN 
    1662                   zmld(:,:) = hmld_kara(:,:) 
    1663                ELSE 
    1664                   CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & 
    1665                      &           ' but ln_kara=.false.' ) 
    1666                ENDIF 
    1667 #else 
    1668                CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & 
    1669                   &           ' but is not defined' ) 
    1670 #endif 
    1671             CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points] 
    1672                !zmld(:,:) = hmld_tref(:,:) 
    1673                CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', & 
    1674                   &           ' but is not available in this version' ) 
    1675             CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points] 
    1676                zmld(:,:) = hmlpt(:,:) 
    1677             END SELECT 
    1678          ENDIF 
    1679  
    1680          ! Propagate through mixed layer 
    1681          DO jj = 1, jpj 
    1682             DO ji = 1, jpi 
    1683                ! 
    1684                jkmax = jpk-1 
    1685                DO jk = jpk-1, 1, -1 
    1686                   IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. & 
    1687                      & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 
    1688                      zmld(ji,jj) = gdepw_n(ji,jj,jk+1) 
    1689                      jkmax = jk 
    1690                   ENDIF 
    1691                END DO 
    1692                ! 
    1693 #if defined key_top 
    1694                DO jk = 2, jkmax 
    1695                   pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:) 
    1696                END DO 
    1697 #endif 
    1698                ! 
    1699             END DO 
    1700          END DO 
    1701  
    1702       ENDIF 
    1703  
    1704       IF ( ln_asmiau ) THEN 
    1705  
    1706          !-------------------------------------------------------------------- 
    1707          ! Incremental Analysis Updating 
    1708          !-------------------------------------------------------------------- 
    1709  
    1710          IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
    1711  
    1712             it = kt - nit000 + 1 
    1713             zincwgt = wgtiau(it)      ! IAU weight for the current time step 
    1714             ! note this is not a tendency so should not be divided by rdt 
    1715  
    1716             IF(lwp) THEN 
    1717                WRITE(numout,*)  
    1718                IF ( ln_pco2inc ) THEN 
    1719                   WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', & 
    1720                      &  kt,' with IAU weight = ', wgtiau(it) 
    1721                ELSE IF ( ln_fco2inc ) THEN 
    1722                   WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', & 
    1723                      &  kt,' with IAU weight = ', wgtiau(it) 
    1724                ENDIF 
    1725                WRITE(numout,*) '~~~~~~~~~~~~' 
    1726             ENDIF 
    1727  
    1728             ! Update the biogeochemical variables 
    1729             ! Add directly to trn and trb, rather than to tra, because tra gets 
    1730             ! reset to zero at the start of trc_stp, called after this routine 
    1731 #if defined key_medusa && defined key_foam_medusa 
    1732             DO jk = 1, jpkm1 
    1733                trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 
    1734                   &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
    1735                trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 
    1736                   &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
    1737             END DO 
    1738 #elif defined key_hadocc 
    1739             DO jk = 1, jpkm1 
    1740                trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 
    1741                   &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
    1742                trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 
    1743                   &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
    1744             END DO 
    1745 #endif 
    1746  
    1747             ! Do not deallocate arrays - needed by asm_bal_wri 
    1748             ! which is called at end of model run 
    1749  
    1750          ENDIF 
    1751  
    1752       ELSEIF ( ln_asmdin ) THEN  
    1753  
    1754          !-------------------------------------------------------------------- 
    1755          ! Direct Initialization 
    1756          !-------------------------------------------------------------------- 
    1757           
    1758          IF ( kt == nitdin_r ) THEN 
    1759  
    1760             neuler = 0                    ! Force Euler forward step 
    1761  
    1762 #if defined key_medusa && defined key_foam_medusa 
    1763             ! Initialize the now fields with the background + increment 
    1764             ! Background currently is what the model is initialised with 
    1765             CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', & 
    1766                &           ' Background state is taken from model rather than background file' ) 
    1767             trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 
    1768                &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) 
    1769             trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 
    1770 #elif defined key_hadocc 
    1771             ! Initialize the now fields with the background + increment 
    1772             ! Background currently is what the model is initialised with 
    1773             CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', & 
    1774                &           ' Background state is taken from model rather than background file' ) 
    1775             trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 
    1776                &                         pco2_balinc(:,:,:,jp_had0:jp_had1) 
    1777             trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 
    1778 #endif 
    1779   
    1780             ! Do not deallocate arrays - needed by asm_bal_wri 
    1781             ! which is called at end of model run 
    1782          ENDIF 
    1783          ! 
    1784       ENDIF 
    1785       ! 
    1786    END SUBROUTINE pco2_asm_inc 
     1194      !! ** Purpose : Apply the biogeochemistry assimilation increments 
     1195      !! 
     1196      !! ** Method  : Call relevant routines in asmbgc 
     1197      !! 
     1198      !! ** Action  : Call relevant routines in asmbgc 
     1199      !! 
     1200      !!---------------------------------------------------------------------- 
     1201      !! 
     1202      INTEGER, INTENT(in   ) :: kt        ! Current time step 
     1203      ! 
     1204      INTEGER                :: icycper   ! Dimension of wgtiau 
     1205      !! 
     1206      !!---------------------------------------------------------------------- 
     1207       
     1208      icycper = SIZE( wgtiau ) 
     1209       
     1210      IF( ln_slchltotinc ) THEN 
     1211         CALL slchltot_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau ) 
     1212      ENDIF 
     1213       
     1214      IF( ln_sfco2inc .OR. ln_spco2inc ) THEN 
     1215         CALL spco2_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, & 
     1216            &                ln_trainc, t_bkginc, s_bkginc ) 
     1217      ENDIF 
     1218 
     1219   END SUBROUTINE bgc_asm_inc 
    17871220    
    17881221   !!====================================================================== 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r9292 r9296  
    6161   USE asminc          ! assimilation increments      
    6262   USE asmbkg          ! writing out state trajectory 
    63    USE asmbal          ! writing out biogeochemical assimilation balancing increments 
     63   USE asmbgc          ! biogeochemical assimilation increments 
    6464   USE diaptr          ! poleward transports           (dia_ptr_init routine) 
    6565   USE diadct          ! sections transports           (dia_dct_init routine) 
     
    161161                IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics 
    162162                IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH 
    163                 IF( ln_logchlinc ) CALL logchl_asm_inc( nit000 - 1 ) 
    164                 IF( ln_fco2inc .OR. ln_pco2inc ) CALL pco2_asm_inc( nit000 - 1 ) 
     163                IF( lk_bgcinc ) CALL bgc_asm_inc( nit000 - 1 )    ! BGC 
    165164             ENDIF 
    166165          ENDIF 
     
    195194      IF( lk_diaobs   )   CALL dia_obs_wri 
    196195      ! 
    197       IF( ( lk_asminc ).AND.( ln_balwri ) ) CALL asm_bal_wri( nitend )  ! Output balancing increments 
     196      IF( ( lk_asminc ).AND.( ln_balwri ) ) CALL asm_bgc_bal_wri( nitend )  ! Output balancing increments 
    198197      ! 
    199198      IF( ln_icebergs )   CALL icb_end( nitend ) 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90

    r9292 r9296  
    253253      ! Passive Tracer Model 
    254254      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    255       IF( lk_asminc .AND. ln_asmiau .AND. ln_logchlinc ) & 
    256          &               CALL logchl_asm_inc( kstp )  ! logchl assimilation 
    257       IF( lk_asminc .AND. ln_asmiau .AND. ( ln_fco2inc .OR. ln_pco2inc ) ) & 
    258          &               CALL pco2_asm_inc( kstp )    ! pCO2/fCO2 assimilation 
     255      IF( lk_asminc .AND. ln_asmiau .AND. lk_bgcinc ) & 
     256         &               CALL bgc_asm_inc( kstp )     ! biogeochemistry assimilation 
    259257                         CALL trc_stp( kstp )         ! time-stepping 
    260258#endif 
Note: See TracChangeset for help on using the changeset viewer.