Changeset 8495


Ignore:
Timestamp:
2017-09-04T18:50:10+02:00 (3 years ago)
Author:
dford
Message:

Merge in changes from dev_r5518_GO6_package_asm_surf_bgc, and adapt to the updated MEDUSA structure.

Location:
branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM
Files:
15 edited
4 copied

Legend:

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

    r8447 r8495  
    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 
    12961300    ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing 
    12971301    nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1] 
     
    13031307    salfixmin = -9999      !  Minimum salinity after applying the increments 
    13041308    nn_divdmp = 0          !  Number of iterations of divergence damping operator 
     1309    mld_choice_bgc = 1     !  MLD criterion to use for biogeochemistry assimilation 
     1310    rn_maxchlinc = -999.0  !  maximum absolute non-log chlorophyll increment from logchl assimilation 
     1311                           !  <= 0 implies no maximum applied (switch turned off) 
     1312                           !   > 0 implies maximum absolute chl increment capped at this value 
    13051313/ 
    13061314!----------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r8400 r8495  
    4949#if defined key_lim3 
    5050   USE ice 
     51#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 
    5168#endif 
    5269   IMPLICIT NONE 
     
    121138!            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
    122139            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) ) 
     173#endif 
    123174            ! 
    124175            CALL iom_close( inum ) 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r8400 r8495  
    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 
    2426   !!---------------------------------------------------------------------- 
    2527   USE wrk_nemo         ! Memory Allocation 
     
    4446#endif 
    4547   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 
    4678 
    4779   IMPLICIT NONE 
     
    5486   PUBLIC   ssh_asm_inc    !: Apply the SSH increment 
    5587   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 
    5690 
    5791#if defined key_asminc 
     
    6195#endif 
    6296   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 
    6398   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment 
    6499   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization 
     
    67102   LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment 
    68103   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 
    69108   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check 
    70109   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
     
    91130   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    92131   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] 
    93160 
    94161   !! * Substitutions 
     
    133200      REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace 
    134201      !! 
    135       NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
     202      NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri,                           & 
    136203         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
     204         &                 ln_logchlinc, ln_logchlbal,                     & 
     205         &                 ln_pco2inc, ln_fco2inc,                         & 
    137206         &                 ln_asmdin, ln_asmiau,                           & 
    138207         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    139208         &                 ln_salfix, salfixmin, nn_divdmp,                & 
    140          &                 ln_seaiceinc, ln_temnofreeze 
     209         &                 ln_seaiceinc, ln_temnofreeze,                   & 
     210         &                 mld_choice_bgc, rn_maxchlinc 
    141211      !!---------------------------------------------------------------------- 
    142212 
     
    163233         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters' 
    164234         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri 
     235         WRITE(numout,*) '      Logical switch for writing out balancing increments      ln_balwri = ', ln_balwri 
    165236         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc 
    166237         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc 
     
    168239         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin 
    169240         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 
    170245         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau 
    171246         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg 
     
    176251         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
    177252         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin 
     253         WRITE(numout,*) '      Choice of MLD for BGC assimilation                  mld_choice_bgc = ', mld_choice_bgc 
     254         WRITE(numout,*) '      Maximum absolute chlorophyll increment (<=0 = off)    rn_maxchlinc = ', rn_maxchlinc 
    178255      ENDIF 
    179256 
     
    223300 
    224301      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 
    225            .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 
    226          & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 
     302         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. & 
     303         &        ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) )) & 
     304         & 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.', & 
    227306         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', & 
    228307         &                ' Inconsistent options') 
     
    233312 
    234313      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 
    235          &                     )  & 
    236          & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 
     314         & .AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) )  & 
     315         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
     316         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc are set to .false. :', & 
    237317         &                ' The assimilation increments are not applied') 
    238318 
     
    254334         &                ' Background time step for Direct Initialization is outside', & 
    255335         &                ' 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 
    256352 
    257353      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further 
     
    350446      ssh_iau(:,:)    = 0.0 
    351447#endif 
    352       IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 
     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 
     464      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 
     465         &  .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN 
    353466 
    354467         !-------------------------------------------------------------------- 
     
    424537         ENDIF 
    425538 
     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 
     559         ENDIF 
     560 
    426561         CALL iom_close( inum ) 
    427562  
     
    535670 
    536671      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 
     693          
     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         ENDIF 
     753 
     754         CALL iom_close( inum ) 
     755          
     756         DO jt = 1, jptra 
     757            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 
     758         END DO 
     759       
     760      ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 
     761 
     762         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) 
     763         tracer_bkg(:,:,:,:) = 0.0 
     764 
     765         CALL iom_open( c_asmbkg, inum ) 
     766          
     767#if defined key_hadocc 
     768         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 
     769         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) 
     770#elif defined key_medusa 
     771         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 
     772         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 
     773#endif 
     774 
     775         CALL iom_close( inum ) 
     776          
     777         DO jt = 1, jptra 
     778            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 
     779         END DO 
     780       
     781      ENDIF 
     782#endif 
    537783      ! 
    538784   END SUBROUTINE asm_inc_init 
     
    11551401 
    11561402   END SUBROUTINE seaice_asm_inc 
     1403 
     1404   SUBROUTINE logchl_asm_inc( kt ) 
     1405      !!---------------------------------------------------------------------- 
     1406      !!                    ***  ROUTINE logchl_asm_inc  *** 
     1407      !!           
     1408      !! ** Purpose : Apply the chlorophyll assimilation increments. 
     1409      !! 
     1410      !! ** Method  : Calculate increments to state variables using nitrogen 
     1411      !!              balancing. 
     1412      !!              Direct initialization or Incremental Analysis Updating. 
     1413      !! 
     1414      !! ** Action  :  
     1415      !!---------------------------------------------------------------------- 
     1416      INTEGER, INTENT(IN) :: kt   ! Current time step 
     1417      ! 
     1418      INTEGER  :: jk              ! Loop counter 
     1419      INTEGER  :: it              ! Index 
     1420      REAL(wp) :: zincwgt         ! IAU weight for current time step 
     1421      REAL(wp) :: zincper         ! IAU interval in seconds 
     1422      !!---------------------------------------------------------------------- 
     1423 
     1424      IF ( kt <= nit000 ) THEN 
     1425 
     1426         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt 
     1427 
     1428#if defined key_medusa && defined key_foam_medusa 
     1429         CALL asm_logchl_bal_medusa( logchl_bkginc, zincper, mld_choice_bgc, & 
     1430            &                        rn_maxchlinc, ln_logchlbal,             & 
     1431            &                        pgrow_avg_bkg, ploss_avg_bkg,           & 
     1432            &                        phyt_avg_bkg, mld_max_bkg,              & 
     1433            &                        tracer_bkg, logchl_balinc ) 
     1434#elif defined key_hadocc 
     1435         CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, & 
     1436            &                        rn_maxchlinc, ln_logchlbal,             & 
     1437            &                        pgrow_avg_bkg, ploss_avg_bkg,           & 
     1438            &                        phyt_avg_bkg, mld_max_bkg,              & 
     1439            &                        chl_bkg, cchl_p_bkg,                    & 
     1440            &                        tracer_bkg, logchl_balinc ) 
     1441#else 
     1442         CALL ctl_stop( 'Attempting to assimilate logchl, ', & 
     1443            &           'but not defined a biogeochemical model' ) 
     1444#endif 
     1445 
     1446      ENDIF 
     1447 
     1448      IF ( ln_asmiau ) THEN 
     1449 
     1450         !-------------------------------------------------------------------- 
     1451         ! Incremental Analysis Updating 
     1452         !-------------------------------------------------------------------- 
     1453 
     1454         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
     1455 
     1456            it = kt - nit000 + 1 
     1457            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
     1458            ! note this is not a tendency so should not be divided by rdt 
     1459 
     1460            IF(lwp) THEN 
     1461               WRITE(numout,*)  
     1462               WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', & 
     1463                  &  kt,' with IAU weight = ', wgtiau(it) 
     1464               WRITE(numout,*) '~~~~~~~~~~~~' 
     1465            ENDIF 
     1466 
     1467            ! Update the biogeochemical variables 
     1468            ! Add directly to trn and trb, rather than to tra, as not a tendency 
     1469#if defined key_medusa && defined key_foam_medusa 
     1470            DO jk = 1, jpkm1 
     1471               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 
     1472                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
     1473               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 
     1474                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
     1475            END DO 
     1476#elif defined key_hadocc 
     1477            DO jk = 1, jpkm1 
     1478               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 
     1479                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
     1480               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 
     1481                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
     1482            END DO 
     1483#endif 
     1484 
     1485            ! Do not deallocate arrays - needed by asm_bal_wri 
     1486            ! which is called at end of model run 
     1487         ENDIF 
     1488 
     1489      ELSEIF ( ln_asmdin ) THEN  
     1490 
     1491         !-------------------------------------------------------------------- 
     1492         ! Direct Initialization 
     1493         !-------------------------------------------------------------------- 
     1494          
     1495         IF ( kt == nitdin_r ) THEN 
     1496 
     1497            neuler = 0                    ! Force Euler forward step 
     1498 
     1499#if defined key_medusa && defined key_foam_medusa 
     1500            ! Initialize the now fields with the background + increment 
     1501            ! Background currently is what the model is initialised with 
     1502            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', & 
     1503               &           ' Background state is taken from model rather than background file' ) 
     1504            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 
     1505               &                         logchl_balinc(:,:,:,jp_msa0:jp_msa1) 
     1506            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 
     1507#elif defined key_hadocc 
     1508            ! Initialize the now fields with the background + increment 
     1509            ! Background currently is what the model is initialised with 
     1510            CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', & 
     1511               &           ' Background state is taken from model rather than background file' ) 
     1512            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 
     1513               &                         logchl_balinc(:,:,:,jp_had0:jp_had1) 
     1514            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 
     1515#endif 
     1516  
     1517            ! Do not deallocate arrays - needed by asm_bal_wri 
     1518            ! which is called at end of model run 
     1519         ENDIF 
     1520         ! 
     1521      ENDIF 
     1522      ! 
     1523   END SUBROUTINE logchl_asm_inc 
     1524 
     1525 
     1526   SUBROUTINE pco2_asm_inc( kt ) 
     1527      !!---------------------------------------------------------------------- 
     1528      !!                    ***  ROUTINE pco2_asm_inc  *** 
     1529      !!           
     1530      !! ** Purpose : Apply the pco2/fco2 assimilation increments. 
     1531      !! 
     1532      !! ** Method  : Calculate increments to state variables using carbon 
     1533      !!              balancing. 
     1534      !!              Direct initialization or Incremental Analysis Updating. 
     1535      !! 
     1536      !! ** Action  :  
     1537      !!---------------------------------------------------------------------- 
     1538      INTEGER, INTENT(IN)                   :: kt           ! Current time step 
     1539      ! 
     1540      INTEGER                               :: ji, jj, jk   ! Loop counters 
     1541      INTEGER                               :: it           ! Index 
     1542      INTEGER                               :: jkmax        ! Index of mixed layer depth 
     1543      REAL(wp)                              :: zincwgt      ! IAU weight for current time step 
     1544      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing 
     1545      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion 
     1546      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background 
     1547      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background 
     1548      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background 
     1549      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background 
     1550      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure 
     1551      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure 
     1552 
     1553      ! Coefficients for fCO2 to pCO2 conversion 
     1554      REAL(wp)                              :: zcoef_fco2_1 = -1636.75 
     1555      REAL(wp)                              :: zcoef_fco2_2 = 12.0408 
     1556      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957 
     1557      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528 
     1558      REAL(wp)                              :: zcoef_fco2_5 = 2.0 
     1559      REAL(wp)                              :: zcoef_fco2_6 = 57.7 
     1560      REAL(wp)                              :: zcoef_fco2_7 = 0.118 
     1561      REAL(wp)                              :: zcoef_fco2_8 = 82.0578 
     1562      !!---------------------------------------------------------------------- 
     1563 
     1564      IF ( kt <= nit000 ) THEN 
     1565 
     1566         !-------------------------------------------------------------------- 
     1567         ! DIC and alkalinity balancing 
     1568         !-------------------------------------------------------------------- 
     1569 
     1570         IF ( ln_fco2inc ) THEN 
     1571#if defined key_medusa && defined key_roam 
     1572            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine 
     1573            patm(1) = 1.0 
     1574            phyd(1) = 0.0 
     1575            DO jj = 1, jpj 
     1576               DO ji = 1, jpi 
     1577                  CALL f2pCO2( pco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) ) 
     1578               END DO 
     1579            END DO 
     1580#else 
     1581            ! If assimilating fCO2, then convert to pCO2 using temperature 
     1582            ! See flux_gas.F90 within HadOCC for details of calculation 
     1583            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) /                                                             & 
     1584               &                    EXP((zcoef_fco2_1                                                            + & 
     1585               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - & 
     1586               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + & 
     1587               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + & 
     1588               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / & 
     1589               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0))) 
     1590#endif 
     1591         ELSE 
     1592            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) 
     1593         ENDIF 
     1594 
     1595         ! Account for physics assimilation if required 
     1596         IF ( ln_trainc ) THEN 
     1597            tem_bkg_temp(:,:) = tsn(:,:,1,1) + t_bkginc(:,:,1) 
     1598            sal_bkg_temp(:,:) = tsn(:,:,1,2) + s_bkginc(:,:,1) 
     1599         ELSE 
     1600            tem_bkg_temp(:,:) = tsn(:,:,1,1) 
     1601            sal_bkg_temp(:,:) = tsn(:,:,1,2) 
     1602         ENDIF 
     1603 
     1604#if defined key_medusa 
     1605         ! Account for logchl balancing if required 
     1606         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN 
     1607            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + logchl_balinc(:,:,1,jpdic) 
     1608            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + logchl_balinc(:,:,1,jpalk) 
     1609         ELSE 
     1610            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) 
     1611            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) 
     1612         ENDIF 
     1613 
     1614         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 
     1615            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        & 
     1616            &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) ) 
     1617 
     1618#elif defined key_hadocc 
     1619         ! Account for logchl balancing if required 
     1620         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN 
     1621            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + logchl_balinc(:,:,1,jp_had_dic) 
     1622            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + logchl_balinc(:,:,1,jp_had_alk) 
     1623         ELSE 
     1624            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) 
     1625            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) 
     1626         ENDIF 
     1627 
     1628         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 
     1629            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        & 
     1630            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) ) 
     1631 
     1632#else 
     1633         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', & 
     1634            &           'but not defined a biogeochemical model' ) 
     1635#endif 
     1636 
     1637         ! Select mixed layer 
     1638         SELECT CASE( mld_choice_bgc ) 
     1639         CASE ( 1 )                   ! Turbocline/mixing depth [W points] 
     1640            zmld(:,:) = hmld(:,:) 
     1641         CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points] 
     1642            zmld(:,:) = hmlp(:,:) 
     1643         CASE ( 3 )                   ! Kara MLD [Interpolated] 
     1644#if defined key_karaml 
     1645            IF ( ln_kara ) THEN 
     1646               zmld(:,:) = hmld_kara(:,:) 
     1647            ELSE 
     1648               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & 
     1649                  &           ' but ln_kara=.false.' ) 
     1650            ENDIF 
     1651#else 
     1652            CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & 
     1653               &           ' but is not defined' ) 
     1654#endif 
     1655         CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points] 
     1656            !zmld(:,:) = hmld_tref(:,:) 
     1657            CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', & 
     1658               &           ' but is not available in this version' ) 
     1659         CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points] 
     1660            zmld(:,:) = hmlpt(:,:) 
     1661         END SELECT 
     1662 
     1663         ! Propagate through mixed layer 
     1664         DO jj = 1, jpj 
     1665            DO ji = 1, jpi 
     1666               ! 
     1667               jkmax = jpk-1 
     1668               DO jk = jpk-1, 1, -1 
     1669                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. & 
     1670                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 
     1671                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1) 
     1672                     jkmax = jk 
     1673                  ENDIF 
     1674               END DO 
     1675               ! 
     1676               DO jk = 2, jkmax 
     1677                  pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:) 
     1678               END DO 
     1679               ! 
     1680            END DO 
     1681         END DO 
     1682 
     1683      ENDIF 
     1684 
     1685      IF ( ln_asmiau ) THEN 
     1686 
     1687         !-------------------------------------------------------------------- 
     1688         ! Incremental Analysis Updating 
     1689         !-------------------------------------------------------------------- 
     1690 
     1691         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
     1692 
     1693            it = kt - nit000 + 1 
     1694            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
     1695            ! note this is not a tendency so should not be divided by rdt 
     1696 
     1697            IF(lwp) THEN 
     1698               WRITE(numout,*)  
     1699               IF ( ln_pco2inc ) THEN 
     1700                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', & 
     1701                     &  kt,' with IAU weight = ', wgtiau(it) 
     1702               ELSE IF ( ln_fco2inc ) THEN 
     1703                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', & 
     1704                     &  kt,' with IAU weight = ', wgtiau(it) 
     1705               ENDIF 
     1706               WRITE(numout,*) '~~~~~~~~~~~~' 
     1707            ENDIF 
     1708 
     1709            ! Update the biogeochemical variables 
     1710            ! Add directly to trn and trb, rather than to tra, as not a tendency 
     1711#if defined key_medusa && defined key_foam_medusa 
     1712            DO jk = 1, jpkm1 
     1713               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 
     1714                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
     1715               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 
     1716                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
     1717            END DO 
     1718#elif defined key_hadocc 
     1719            DO jk = 1, jpkm1 
     1720               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 
     1721                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
     1722               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 
     1723                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
     1724            END DO 
     1725#endif 
     1726 
     1727            ! Do not deallocate arrays - needed by asm_bal_wri 
     1728            ! which is called at end of model run 
     1729 
     1730         ENDIF 
     1731 
     1732      ELSEIF ( ln_asmdin ) THEN  
     1733 
     1734         !-------------------------------------------------------------------- 
     1735         ! Direct Initialization 
     1736         !-------------------------------------------------------------------- 
     1737          
     1738         IF ( kt == nitdin_r ) THEN 
     1739 
     1740            neuler = 0                    ! Force Euler forward step 
     1741 
     1742#if defined key_medusa && defined key_foam_medusa 
     1743            ! Initialize the now fields with the background + increment 
     1744            ! Background currently is what the model is initialised with 
     1745            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', & 
     1746               &           ' Background state is taken from model rather than background file' ) 
     1747            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 
     1748               &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) 
     1749            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 
     1750#elif defined key_hadocc 
     1751            ! Initialize the now fields with the background + increment 
     1752            ! Background currently is what the model is initialised with 
     1753            CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', & 
     1754               &           ' Background state is taken from model rather than background file' ) 
     1755            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 
     1756               &                         pco2_balinc(:,:,:,jp_had0:jp_had1) 
     1757            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 
     1758#endif 
     1759  
     1760            ! Do not deallocate arrays - needed by asm_bal_wri 
     1761            ! which is called at end of model run 
     1762         ENDIF 
     1763         ! 
     1764      ENDIF 
     1765      ! 
     1766   END SUBROUTINE pco2_asm_inc 
    11571767    
    11581768   !!====================================================================== 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90

    r6486 r8495  
    2020      & c_asmtrj = 'assim_trj',                  & !: Filename for storing the  
    2121                                                   !: reference trajectory 
    22       & c_asminc = 'assim_background_increments'  !: Filename for storing the  
     22      & c_asminc = 'assim_background_increments', & !: Filename for storing the  
    2323                                                   !: increments to the background 
    2424                                                   !: state 
     25      & c_asmbal = 'assim.balincs'                 !: Filename for storing the  
     26                                                   !: balancing increments calculated 
     27                                                   !: for biogeochemistry 
    2528 
    2629   INTEGER, PUBLIC :: nitbkg_r      !: Background time step referenced to nit000 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r7179 r8495  
    3838   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: tidal mixing flag 
    3939 
    40    !                       !!* Namelist  namzdf_tmx : tidal mixing * 
    41    REAL(wp) ::  rn_htmx     ! vertical decay scale for turbulence (meters) 
    42    REAL(wp) ::  rn_n2min    ! threshold of the Brunt-Vaisala frequency (s-1) 
    43    REAL(wp) ::  rn_tfe      ! tidal dissipation efficiency (St Laurent et al. 2002) 
    44    REAL(wp) ::  rn_me       ! mixing efficiency (Osborn 1980) 
    45    LOGICAL  ::  ln_tmx_itf  ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization 
    46    REAL(wp) ::  rn_tfe_itf  ! ITF tidal dissipation efficiency (St Laurent et al. 2002) 
    47  
    48    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   en_tmx     ! energy available for tidal mixing (W/m2) 
    49    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)  ::   mask_itf   ! mask to use over Indonesian area 
    50    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   az_tmx     ! coefficient used to evaluate the tidal induced Kz 
     40   !                              !!* Namelist  namzdf_tmx : tidal mixing * 
     41   REAL(wp)        ::  rn_htmx     ! vertical decay scale for turbulence (meters) 
     42   REAL(wp)        ::  rn_n2min    ! threshold of the Brunt-Vaisala frequency (s-1) 
     43   REAL(wp)        ::  rn_tfe      ! tidal dissipation efficiency (St Laurent et al. 2002) 
     44   REAL(wp)        ::  rn_me       ! mixing efficiency (Osborn 1980) 
     45   LOGICAL, PUBLIC ::  ln_tmx_itf  ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization 
     46   REAL(wp)        ::  rn_tfe_itf  ! ITF tidal dissipation efficiency (St Laurent et al. 2002) 
     47 
     48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   en_tmx     ! energy available for tidal mixing (W/m2) 
     49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mask_itf   ! mask to use over Indonesian area 
     50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)       ::   az_tmx     ! coefficient used to evaluate the tidal induced Kz 
    5151 
    5252   !! * Substitutions 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r8400 r8495  
    6161   USE asminc          ! assimilation increments      
    6262   USE asmbkg          ! writing out state trajectory 
     63   USE asmbal          ! writing out assimilation balancing increments 
    6364   USE diaptr          ! poleward transports           (dia_ptr_init routine) 
    6465   USE diadct          ! sections transports           (dia_dct_init routine) 
     
    160161                IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics 
    161162                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 ) 
    162165             ENDIF 
    163166          ENDIF 
     
    191194 
    192195      IF( lk_diaobs   )   CALL dia_obs_wri 
     196      ! 
     197      IF( ( lk_asminc ).AND.( ln_balwri ) ) CALL asm_bal_wri( nitend )  ! Output balancing increments 
    193198      ! 
    194199      IF( ln_icebergs )   CALL icb_end( nitend ) 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8400 r8495  
    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 
    255259                         CALL trc_stp( kstp )         ! time-stepping 
    256260#endif 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90

    r8442 r8495  
    7676                                   zn_dms_chd, zn_dms_chn, zn_dms_din,    & 
    7777                                   zn_dms_mld, zn_dms_qsr,                & 
     78# if defined key_foam_medusa 
     79                                   f2_pco2w, f2_fco2w,                    & 
     80# endif 
    7881                                   xnln, xnld  
    7982      USE trc,               ONLY: med_diag 
     
    8689#  else 
    8790      USE trcco2_medusa,     ONLY: trc_co2_medusa 
     91#   if defined key_foam_medusa 
     92      USE mocsy_mainmod,     ONLY: p2fCO2 
     93#   endif 
    8894#  endif 
    8995      USE trcdms_medusa,     ONLY: trc_dms_medusa 
     
    234240                     iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt 
    235241               endif 
     242#    if defined key_foam_medusa 
     243               !! DAF (Aug 2017): calculate fCO2 for observation operator 
     244               CALL p2fCO2( f_pco2w, ztmp, f_pp0, 0.0, 1, f_fco2w ) 
     245#    endif 
    236246            ENDIF 
    237247         ENDDO 
     
    408418                              CO2flux_conv 
    409419               !! ENDIF 
     420#   if defined key_foam_medusa 
     421               !! DAF (Aug 2017): Save pCO2 and fCO2 for observation operator 
     422               f2_pco2w(ji,jj) = f_pco2w(ji,jj) 
     423               f2_fco2w(ji,jj) = f_pco2w(ji,jj) 
     424#   endif 
    410425               IF ( lk_iomput ) THEN 
    411426                  IF( med_diag%ATM_PCO2%dgsave ) THEN 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90

    r8442 r8495  
    4545                                   f3_co3, f3_h2co3, f3_hco3,           & 
    4646                                   f3_omarg, f3_omcal, f3_pH,           & 
     47# if defined key_foam_medusa 
     48                                   mld_max, pgrow_avg,                  & 
     49                                   ploss_avg, phyt_avg,                 & 
     50# endif 
    4751                                   za_sed_c, za_sed_ca, za_sed_fe,      & 
    4852                                   za_sed_n, za_sed_si,                 & 
     
    5357      USE trc,               ONLY: med_diag, nittrc000  
    5458      USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd 
     59# if defined key_foam_medusa 
     60      USE zdfmxl,            ONLY: hmld 
     61# endif 
    5562  
    5663      !! time (integer timestep) 
     
    218225      endif 
    219226 
     227# if defined key_foam_medusa 
     228!!---------------------------------------------------------------------- 
     229!! Diagnostics required for ocean colour assimilation: 
     230!! Mixed layer average phytoplankton growth, loss and concentration 
     231!! Maximum mixed layer depth 
     232!!---------------------------------------------------------------------- 
     233!! 
     234      DO jj = 2,jpjm1 
     235         DO ji = 2,jpim1 
     236            pgrow_avg(ji,jj) = pgrow_avg(ji,jj) / hmld(ji,jj) 
     237            ploss_avg(ji,jj) = ploss_avg(ji,jj) / hmld(ji,jj) 
     238            phyt_avg(ji,jj)  = phyt_avg(ji,jj)  / hmld(ji,jj) 
     239            IF ( hmld(ji,jj) .GT. mld_max(ji,jj) ) THEN 
     240               mld_max(ji,jj) = hmld(ji,jj) 
     241            ENDIF 
     242         END DO 
     243      END DO 
     244# endif 
     245 
    220246#  if defined key_debug_medusa 
    221247         !! AXY (12/07/17) 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90

    r8442 r8495  
    3434      USE bio_medusa_mod 
    3535      USE par_oce,           ONLY: jpi, jpj, jpk 
     36# if defined key_foam_medusa 
     37      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max 
     38# else 
    3639      USE sms_medusa,        ONLY: jdms 
     40# endif 
    3741      USE trc,               ONLY: ln_diatrc, med_diag, nittrc000  
    3842      USE in_out_manager,    ONLY: lwp 
     
    164168      fslownflux(:,:) = 0.0  
    165169      fslowcflux(:,:) = 0.0  
     170      !! 
     171# if defined key_foam_medusa 
     172      pgrow_avg(:,:) = 0.0 
     173      ploss_avg(:,:) = 0.0 
     174      phyt_avg(:,:)  = 0.0 
     175      IF( kt == nittrc000 ) THEN 
     176         mld_max(:,:) = 0.0 
     177      ENDIF 
     178# endif 
    166179      !! 
    167180      !! allocate and initiate 2D diag 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90

    r8441 r8495  
    3636                                   fdpn, fdpn2, fdzme, fdzme2,             & 
    3737                                   fdzmi, fdzmi2, fsdiss, fsin,            & 
     38# if defined key_foam_medusa 
     39                                   fdep1, fprn, fprd,                      & 
     40                                   fgmepd, fgmepn, fgmipn,                 & 
     41# endif 
    3842                                   zphd, zphn, zpds, zzme, zzmi 
     43# if defined key_foam_medusa 
     44      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask 
     45      USE par_kind,          ONLY: wp 
     46# else 
    3947      USE dom_oce,           ONLY: tmask 
     48# endif 
    4049      USE par_oce,           ONLY: jpim1, jpjm1 
    4150      USE phytoplankton_mod, ONLY: phytoplankton 
    4251      USE sms_medusa,        ONLY: jmpd, jmpn, jmzme, jmzmi,               & 
     52# if defined key_foam_medusa 
     53                                   pgrow_avg, ploss_avg, phyt_avg,         & 
     54# endif 
    4355                                   xkphd, xkphn, xkzme, xkzmi,             & 
    4456                                   xmetapd, xmetapn, xmetazme, xmetazmi,   & 
    4557                                   xmpd, xmpn, xmzme, xmzmi, xsdiss 
     58# if defined key_foam_medusa 
     59      USE zdfmxl,            ONLY: hmld 
     60# endif 
    4661      USE zooplankton_mod,   ONLY: zooplankton 
     62 
     63# if defined key_foam_medusa 
     64   !!* Substitution 
     65#  include "domzgr_substitute.h90" 
     66# endif 
    4767 
    4868      !! Level 
     
    5070 
    5171      INTEGER :: ji, jj 
     72 
     73# if defined key_foam_medusa 
     74      REAL(wp) :: fq0 
     75# endif 
    5276 
    5377      !!------------------------------------------------------------------- 
     
    188212      ENDDO 
    189213 
     214# if defined key_foam_medusa 
     215      !! Mixed layer averages for ocean colour assimilation 
     216      !! 
     217      if (fdep1(ji,jj).le.hmld(ji,jj)) then 
     218         !! this level is entirely in the mixed layer 
     219         fq0 = 1.0 
     220      elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then 
     221         !! this level is entirely below the mixed layer 
     222         fq0 = 0.0 
     223      else 
     224         !! this level straddles the mixed layer 
     225         fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk) 
     226      endif 
     227      !! 
     228      pgrow_avg(ji,jj) = pgrow_avg(ji,jj) +                                  & 
     229                         (((fprn(ji,jj) * zphn(ji,jj)) +                     & 
     230                           (fprd(ji,jj) * zphd(ji,jj))  ) *                  & 
     231                          fse3t(ji,jj,jk) * fq0) 
     232      ploss_avg(ji,jj) = ploss_avg(ji,jj) +  & 
     233                         ((fgmepd(ji,jj) + fdpd(ji,jj) + fdpd2(ji,jj) +      & 
     234                           fgmepn(ji,jj) + fdpn(ji,jj) + fdpn2(ji,jj) +      & 
     235                           fgmipn(ji,jj)                               ) *   & 
     236                          fse3t(ji,jj,jk) * fq0) 
     237      phyt_avg(ji,jj)  = phyt_avg(ji,jj)  +                                  & 
     238                         ((zphn(ji,jj) + zphd(ji,jj)) *                      & 
     239                          fse3t(ji,jj,jk) * fq0) 
     240# endif 
     241 
    190242   END SUBROUTINE plankton 
    191243 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90

    r8132 r8495  
    205205   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_ccd_arg  !: 2D aragonite CCD depth 
    206206!! 
     207#if defined key_foam_medusa 
     208!! 2D fields of pCO2 and fCO2 for observation operator 
     209   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_pco2w    !: 2D pCO2 
     210   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_fco2w    !: 2D fCO2 
     211!! 
     212#endif 
    207213!! 2D fields of organic and inorganic material sedimented on the seafloor 
    208214   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_sed_n    !: 2D organic nitrogen   (before) 
     
    384390      & 926.4960, 935.7040 / 
    385391# endif 
     392# if defined key_foam_medusa 
     393   REAL(wp) ::   xobs_xco2a   !: Observed atmospheric xCO2, read in 
     394# endif 
    386395#endif 
    387396 
     
    434443   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: cmask       !: ??? 
    435444 
     445#if defined key_foam_medusa 
     446!!---------------------------------------------------------------------- 
     447!! Parameters required for ocean colour assimilation 
     448!!---------------------------------------------------------------------- 
     449!! 
     450   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pgrow_avg  !: Mixed layer average phytoplankton growth 
     451   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ploss_avg  !: Mixed layer average phytoplankton loss 
     452   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: phyt_avg   !: Mixed layer average phytoplankton 
     453   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_max    !: Maximum mixed layer depth 
     454!! 
     455#endif 
     456 
    436457   !!---------------------------------------------------------------------- 
    437458   !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     
    446467      !!---------------------------------------------------------------------- 
    447468      USE lib_mpp , ONLY: ctl_warn 
    448       INTEGER ::   ierr(8)        ! Local variables 
     469      INTEGER ::   ierr(9)        ! Local variables 
    449470      !!---------------------------------------------------------------------- 
    450471      ierr(:) = 0 
     
    456477      !* 2D and 3D fields of carbonate system parameters 
    457478      ALLOCATE( f2_ccd_cal(jpi,jpj)  , f2_ccd_arg(jpi,jpj)  ,       & 
     479#  if defined key_foam_medusa 
     480                f2_pco2w(jpi,jpj)    , f2_fco2w(jpi,jpj)    ,       & 
     481#  endif 
    458482         &      f3_pH(jpi,jpj,jpk)   , f3_h2co3(jpi,jpj,jpk),       & 
    459483         &      f3_hco3(jpi,jpj,jpk) , f3_co3(jpi,jpj,jpk)  ,       & 
     
    504528         &      ffln(jpi,jpj,jpk)    , fflf(jpi,jpj,jpk)    ,       & 
    505529         &      ffls(jpi,jpj,jpk)    , cmask(jpi,jpj)       ,    STAT=ierr(8) )  
     530# if defined key_foam_medusa 
     531      ALLOCATE( pgrow_avg(jpi,jpj)   , ploss_avg(jpi,jpj)   ,       & 
     532         &      phyt_avg(jpi,jpj)    , mld_max(jpi,jpj)     ,    STAT=ierr(9) ) 
     533# endif 
    506534#endif 
    507535      ! 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r8442 r8495  
    9898      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm 
    9999      USE sbc_oce,                    ONLY: lk_oasis 
     100# if defined key_foam_medusa 
     101      USE sms_medusa,                 ONLY: hist_pco2, xobs_xco2a,          & 
     102                                            pgrow_avg, ploss_avg,           & 
     103                                            phyt_avg, mld_max 
     104# else 
    100105      USE sms_medusa,                 ONLY: hist_pco2 
     106# endif 
    101107      USE trc,                        ONLY: ln_rsttr, nittrc000, trn 
    102108      USE bio_medusa_init_mod,        ONLY: bio_medusa_init 
     
    309315         f_xco2a(:,:) = fq4 
    310316      endif 
     317#  if defined key_foam_medusa 
     318      IF ( xobs_xco2a > 0.0 ) THEN 
     319         IF(lwp) WRITE(numout,*) ' using observed atm pCO2 = ', xobs_xco2a 
     320         f_xco2a(:,:) = xobs_xco2a 
     321      ELSE 
     322         IF(lwp) WRITE(numout,*) ' xobs_xco2a <= 0 so using default atm pCO2' 
     323      ENDIF 
     324#  endif 
    311325#  if defined key_axy_pi_co2 
    312326      !! OCMIP pre-industrial pCO2 
     
    348362      !!          x * 30d + 1*rdt(i.e: mod = rdt)    
    349363      !!          ++ need to pass carb-chem output var through restarts 
     364#if defined key_foam_medusa 
     365      !! DAF (Aug 2017): For FOAM we want to run daily 
     366      If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
     367           (mod(kt*rdt,86400.) == rdt) ) THEN 
     368#else 
    350369      If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
    351370           ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN 
    352          !!--------------------------------------------------------------- 
     371#endif 
     372         !!---------------------------------------------------------------------- 
    353373         !! Calculate the carbonate chemistry for the whole ocean on the first 
    354374         !! simulation timestep and every month subsequently; the resulting 3D 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90

    r8442 r8495  
    8888#if defined key_roam 
    8989      NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        & 
     90# if defined key_foam_medusa 
     91      &    xobs_xco2a,                                        & 
     92# endif 
    9093      &    xthetarem,xo2min  
    9194#endif 
     
    474477         &   ' key_debug_medusa                                                       = INACTIVE' 
    475478# endif 
     479#if defined key_foam_medusa 
     480         WRITE(numout,*)     & 
     481         &   ' key_foam_medusa                                                        = ACTIVE' 
     482#else 
     483         WRITE(numout,*)     & 
     484         &   ' key_foam_medusa                                                        = INACTIVE' 
     485#endif 
    476486         WRITE(numout,*) ' ' 
    477487 
     
    9901000      xthetarem = 0. 
    9911001      xo2min    = 0. 
     1002# if defined key_foam_medusa 
     1003      xobs_xco2a = 0. 
     1004# endif 
    9921005 
    9931006      !READ(numnatm,natroam) 
     
    10091022!!       xthetarem :  oxygen consumption by carbon remineralisation 
    10101023!!       xo2min    :  oxygen minimum concentration 
     1024# if defined key_foam_medusa 
     1025!!       xobs_xco2a : observed atmospheric xCO2 
     1026# endif 
    10111027 
    10121028      IF(lwp) THEN 
     
    10261042          WRITE(numout,*)     & 
    10271043          &   ' oxygen minimum concentration                               xo2min      = ', xo2min 
     1044# if defined key_foam_medusa 
     1045          WRITE(numout,*)     & 
     1046          &   ' observed atmospheric xCO2                                  xobs_xco2a  = ', xobs_xco2a 
     1047# endif 
    10281048       ENDIF 
    10291049 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/TOP_SRC/trcrst.F90

    r8427 r8495  
    4343   USE sbc_oce, ONLY: lk_oasis  
    4444   USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable 
     45#if defined key_foam_medusa 
     46   USE obs_const, ONLY: obfillflt  ! Observation operator fill value 
     47#endif 
    4548 
    4649   IMPLICIT NONE 
     
    329332         IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...'  
    330333      ENDIF 
     334      ! 
     335#  if defined key_foam_medusa 
     336      !! 2D fields of pCO2 and fCO2 for observation operator on first timestep 
     337      IF( iom_varid( numrtr, 'PCO2W', ldstop = .FALSE. ) > 0 ) THEN 
     338         IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 present - reading in ...' 
     339         CALL iom_get( numrtr, jpdom_autoglo, 'PCO2W',  f2_pco2w(:,:)  ) 
     340         CALL iom_get( numrtr, jpdom_autoglo, 'FCO2W',  f2_fco2w(:,:)  ) 
     341      ELSE 
     342         IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 absent - setting to fill ...' 
     343         f2_pco2w(:,:) = obfillflt * tmask(:,:,1) 
     344         f2_fco2w(:,:) = obfillflt * tmask(:,:,1) 
     345      ENDIF 
     346#  endif 
    331347# endif 
    332  
     348# if defined key_foam_medusa 
     349      !! Fields for ocean colour assimilation on first timestep 
     350      IF( iom_varid( numrtr, 'pgrow_avg', ldstop = .FALSE. ) > 0 ) THEN 
     351         IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg present - reading in ...' 
     352         CALL iom_get( numrtr, jpdom_autoglo, 'pgrow_avg',  pgrow_avg(:,:)  ) 
     353         CALL iom_get( numrtr, jpdom_autoglo, 'ploss_avg',  ploss_avg(:,:)  ) 
     354         CALL iom_get( numrtr, jpdom_autoglo, 'phyt_avg',   phyt_avg(:,:)   ) 
     355         CALL iom_get( numrtr, jpdom_autoglo, 'mld_max',    mld_max(:,:)    ) 
     356      ELSE 
     357         IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg absent - setting to zero ...' 
     358         pgrow_avg(:,:) = 0.0 
     359         ploss_avg(:,:) = 0.0 
     360         phyt_avg(:,:)  = 0.0 
     361         mld_max(:,:)   = 0.0 
     362      ENDIF 
     363# endif 
    333364 
    334365#endif 
     
    498529      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG') 
    499530      !! 
     531#  if defined key_foam_medusa 
     532      !! Fields for observation operator on first timestep 
     533      IF(lwp) WRITE(numout,*) ' MEDUSA OBS fields - writing out ...' 
     534      CALL iom_rstput( kt, nitrst, numrtw, 'PCO2W', f2_pco2w(:,:)  ) 
     535      CALL iom_rstput( kt, nitrst, numrtw, 'FCO2W', f2_fco2w(:,:)  ) 
     536#  endif 
     537# endif 
     538# if defined key_foam_medusa 
     539      !! Fields for assimilation on first timestep 
     540      IF(lwp) WRITE(numout,*) ' MEDUSA ASM fields - writing out ...' 
     541      CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg', pgrow_avg(:,:) ) 
     542      CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg', ploss_avg(:,:) ) 
     543      CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg',  phyt_avg(:,:)  ) 
     544      CALL iom_rstput( kt, nitrst, numrtw, 'mld_max',   mld_max(:,:)   ) 
    500545# endif 
    501546!! 
Note: See TracChangeset for help on using the changeset viewer.