Changeset 9292


Ignore:
Timestamp:
2018-01-30T19:41:58+01:00 (2 years ago)
Author:
dford
Message:

Merge in changes from dev_r5518_GO6_package_asm_surf_bgc_v2.

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

Legend:

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

    r8447 r9292  
    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_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r8400 r9292  
    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_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r8400 r9292  
    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            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 
    537789      ! 
    538790   END SUBROUTINE asm_inc_init 
     
    11551407 
    11561408   END SUBROUTINE seaice_asm_inc 
     1409 
     1410   SUBROUTINE logchl_asm_inc( kt ) 
     1411      !!---------------------------------------------------------------------- 
     1412      !!                    ***  ROUTINE logchl_asm_inc  *** 
     1413      !!           
     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 
    11571787    
    11581788   !!====================================================================== 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90

    r6486 r9292  
    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_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r7179 r9292  
    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_3d_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r8400 r9292  
    6161   USE asminc          ! assimilation increments      
    6262   USE asmbkg          ! writing out state trajectory 
     63   USE asmbal          ! writing out biogeochemical 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_3d_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8400 r9292  
    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_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90

    r9257 r9292  
    7777                                   zn_dms_chd, zn_dms_chn, zn_dms_din,    & 
    7878                                   zn_dms_mld, zn_dms_qsr,                & 
     79# if defined key_foam_medusa 
     80                                   f2_pco2w, f2_fco2w,                    & 
     81# endif 
    7982                                   xnln, xnld  
    8083      USE trc,               ONLY: med_diag 
     
    8790#  else 
    8891      USE trcco2_medusa,     ONLY: trc_co2_medusa 
     92#   if defined key_foam_medusa 
     93      USE mocsy_mainmod,     ONLY: p2fCO2 
     94#   endif 
    8995#  endif 
    9096      USE trcdms_medusa,     ONLY: trc_dms_medusa 
     
    330336                     iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt 
    331337               endif 
     338#    if defined key_foam_medusa 
     339               !! DAF (Aug 2017): calculate fCO2 for observation operator 
     340               CALL p2fCO2( f_pco2w, ztmp, f_pp0, 0.0, 1, f_fco2w ) 
     341#    endif 
    332342            ENDIF 
    333343         ENDDO 
     
    509519                              CO2flux_conv 
    510520               !! ENDIF 
     521#   if defined key_foam_medusa 
     522               !! DAF (Aug 2017): Save pCO2 and fCO2 for observation operator 
     523               f2_pco2w(ji,jj) = f_pco2w(ji,jj) 
     524               f2_fco2w(ji,jj) = f_pco2w(ji,jj) 
     525#   endif 
    511526               IF ( lk_iomput ) THEN 
    512527                  IF( med_diag%ATM_PCO2%dgsave ) THEN 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90

    r9114 r9292  
    4848                                   f3_co3, f3_h2co3, f3_hco3,           & 
    4949                                   f3_omarg, f3_omcal, f3_pH,           & 
     50# if defined key_foam_medusa 
     51                                   mld_max, pgrow_avg,                  & 
     52                                   ploss_avg, phyt_avg,                 & 
     53# endif 
    5054                                   za_sed_c, za_sed_ca, za_sed_fe,      & 
    5155                                   za_sed_n, za_sed_si,                 & 
     
    5761      USE trc,               ONLY: med_diag, nittrc000, trn  
    5862      USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd 
     63# if defined key_foam_medusa 
     64      USE zdfmxl,            ONLY: hmld 
     65# endif 
    5966  
    6067      !! time (integer timestep) 
     
    115122# endif       
    116123 
     124# if defined key_foam_medusa 
     125!!---------------------------------------------------------------------- 
     126!! Diagnostics required for ocean colour assimilation: 
     127!! Mixed layer average phytoplankton growth, loss and concentration 
     128!! Maximum mixed layer depth 
     129!!---------------------------------------------------------------------- 
     130!! 
     131      DO jj = 2,jpjm1 
     132         DO ji = 2,jpim1 
     133            IF ( hmld(ji,jj) .GT. 0.0 ) THEN 
     134               pgrow_avg(ji,jj) = pgrow_avg(ji,jj) / hmld(ji,jj) 
     135               ploss_avg(ji,jj) = ploss_avg(ji,jj) / hmld(ji,jj) 
     136               phyt_avg(ji,jj)  = phyt_avg(ji,jj)  / hmld(ji,jj) 
     137               IF ( hmld(ji,jj) .GT. mld_max(ji,jj) ) THEN 
     138                  mld_max(ji,jj) = hmld(ji,jj) 
     139               ENDIF 
     140            ENDIF 
     141         END DO 
     142      END DO 
     143# endif 
     144 
    117145#  if defined key_debug_medusa 
    118146         !! AXY (12/07/17) 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90

    r9257 r9292  
    3535      USE bio_medusa_mod 
    3636      USE par_oce,           ONLY: jpi, jpj, jpk 
     37# if defined key_foam_medusa 
     38      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max 
     39# else 
    3740      USE sms_medusa,        ONLY: jdms 
     41# endif 
    3842      USE trc,               ONLY: ln_diatrc, med_diag, nittrc000  
    3943      USE in_out_manager,    ONLY: lwp 
     
    192196      fslowsinkc(:,:) = 0.0 
    193197# endif       
     198      !! 
     199# if defined key_foam_medusa 
     200      pgrow_avg(:,:) = 0.0 
     201      ploss_avg(:,:) = 0.0 
     202      phyt_avg(:,:)  = 0.0 
     203      IF( kt == nittrc000 ) THEN 
     204         mld_max(:,:) = 0.0 
     205      ENDIF 
     206# endif 
    194207      !! 
    195208      !! allocate and initiate 2D diag 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90

    r8441 r9292  
    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      DO jj = 2,jpjm1 
     218         DO ji = 2,jpim1 
     219            IF (tmask(ji,jj,jk) == 1) THEN 
     220               if (fdep1(ji,jj).le.hmld(ji,jj)) then 
     221                  !! this level is entirely in the mixed layer 
     222                  fq0 = 1.0 
     223               elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then 
     224                  !! this level is entirely below the mixed layer 
     225                  fq0 = 0.0 
     226               else 
     227                  !! this level straddles the mixed layer 
     228                  fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk) 
     229               endif 
     230               !! 
     231               pgrow_avg(ji,jj) = pgrow_avg(ji,jj) +                   & 
     232                                  (((fprn(ji,jj) * zphn(ji,jj)) +      & 
     233                                    (fprd(ji,jj) * zphd(ji,jj))  ) *   & 
     234                                   fse3t(ji,jj,jk) * fq0) 
     235               ploss_avg(ji,jj) = ploss_avg(ji,jj) +                   & 
     236                                  ((fgmepd(ji,jj) + fdpd(ji,jj) +      & 
     237                                    fdpd2(ji,jj)                +      & 
     238                                    fgmepn(ji,jj) + fdpn(ji,jj) +      & 
     239                                    fdpn2(ji,jj)  + fgmipn(ji,jj) ) *  & 
     240                                   fse3t(ji,jj,jk) * fq0) 
     241               phyt_avg(ji,jj)  = phyt_avg(ji,jj)  +                   & 
     242                                  ((zphn(ji,jj) + zphd(ji,jj)) *       & 
     243                                   fse3t(ji,jj,jk) * fq0) 
     244            ENDIF 
     245         ENDDO 
     246      ENDDO 
     247# endif 
     248 
    190249   END SUBROUTINE plankton 
    191250 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90

    r9258 r9292  
    212212   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_ccd_arg  !: 2D aragonite CCD depth 
    213213!! 
     214#if defined key_foam_medusa 
     215!! 2D fields of pCO2 and fCO2 for observation operator 
     216   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_pco2w    !: 2D pCO2 
     217   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_fco2w    !: 2D fCO2 
     218!! 
     219#endif 
    214220!! 2D fields of organic and inorganic material sedimented on the seafloor 
    215221   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_sed_n    !: 2D organic nitrogen   (before) 
     
    391397      & 926.4960, 935.7040 / 
    392398# endif 
     399# if defined key_foam_medusa 
     400   REAL(wp) ::   xobs_xco2a   !: Observed atmospheric xCO2, read in 
     401# endif 
    393402#endif 
    394403 
     
    441450   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: cmask       !: ??? 
    442451 
     452#if defined key_foam_medusa 
     453!!---------------------------------------------------------------------- 
     454!! Parameters required for ocean colour assimilation 
     455!!---------------------------------------------------------------------- 
     456!! 
     457   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pgrow_avg  !: Mixed layer average phytoplankton growth 
     458   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ploss_avg  !: Mixed layer average phytoplankton loss 
     459   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: phyt_avg   !: Mixed layer average phytoplankton 
     460   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_max    !: Maximum mixed layer depth 
     461!! 
     462#endif 
     463 
    443464   !!---------------------------------------------------------------------- 
    444465   !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     
    453474      !!---------------------------------------------------------------------- 
    454475      USE lib_mpp , ONLY: ctl_warn 
    455       INTEGER ::   ierr(8)        ! Local variables 
     476      INTEGER ::   ierr(9)        ! Local variables 
    456477      !!---------------------------------------------------------------------- 
    457478      ierr(:) = 0 
     
    463484      !* 2D and 3D fields of carbonate system parameters 
    464485      ALLOCATE( f2_ccd_cal(jpi,jpj)  , f2_ccd_arg(jpi,jpj)  ,       & 
     486#  if defined key_foam_medusa 
     487                f2_pco2w(jpi,jpj)    , f2_fco2w(jpi,jpj)    ,       & 
     488#  endif 
    465489         &      f3_pH(jpi,jpj,jpk)   , f3_h2co3(jpi,jpj,jpk),       & 
    466490         &      f3_hco3(jpi,jpj,jpk) , f3_co3(jpi,jpj,jpk)  ,       & 
     
    511535         &      ffln(jpi,jpj,jpk)    , fflf(jpi,jpj,jpk)    ,       & 
    512536         &      ffls(jpi,jpj,jpk)    , cmask(jpi,jpj)       ,    STAT=ierr(8) )  
     537# if defined key_foam_medusa 
     538      ALLOCATE( pgrow_avg(jpi,jpj)   , ploss_avg(jpi,jpj)   ,       & 
     539         &      phyt_avg(jpi,jpj)    , mld_max(jpi,jpj)     ,    STAT=ierr(9) ) 
     540# endif 
    513541#endif 
    514542      ! 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r9257 r9292  
    102102      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm 
    103103      USE sbc_oce,                    ONLY: lk_oasis 
     104# if defined key_foam_medusa 
     105      USE sms_medusa,                 ONLY: hist_pco2, xobs_xco2a,          & 
     106                                            pgrow_avg, ploss_avg,           & 
     107                                            phyt_avg, mld_max 
     108# else 
    104109      USE sms_medusa,                 ONLY: hist_pco2 
     110# endif 
    105111      USE trc,                        ONLY: ln_rsttr, nittrc000, trn 
    106112      USE bio_medusa_init_mod,        ONLY: bio_medusa_init 
     
    319325         f_xco2a(:,:) = fq4 
    320326      endif 
     327#  if defined key_foam_medusa 
     328      IF ( xobs_xco2a > 0.0 ) THEN 
     329         IF(lwp) WRITE(numout,*) ' using observed atm pCO2 = ', xobs_xco2a 
     330         f_xco2a(:,:) = xobs_xco2a 
     331      ELSE 
     332         IF(lwp) WRITE(numout,*) ' xobs_xco2a <= 0 so using default atm pCO2' 
     333      ENDIF 
     334#  endif 
    321335#  if defined key_axy_pi_co2 
    322336      !! OCMIP pre-industrial pCO2 
     
    358372      !!          x * 30d + 1*rdt(i.e: mod = rdt)    
    359373      !!          ++ need to pass carb-chem output var through restarts 
     374#if defined key_foam_medusa 
     375      !! DAF (Aug 2017): For FOAM we want to run daily 
     376      If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
     377           (mod(kt*rdt,86400.) == rdt) ) THEN 
     378#else 
    360379      If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
    361380           ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN 
    362          !!--------------------------------------------------------------- 
     381#endif 
     382         !!---------------------------------------------------------------------- 
    363383         !! Calculate the carbonate chemistry for the whole ocean on the first 
    364384         !! simulation timestep and every month subsequently; the resulting 3D 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90

    r9258 r9292  
    9090#if defined key_roam 
    9191      NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        & 
     92# if defined key_foam_medusa 
     93      &    xobs_xco2a,                                        & 
     94# endif 
    9295      &    xthetarem,xo2min  
    9396#endif 
     
    490493         &   ' key_debug_medusa                                                       = INACTIVE' 
    491494# endif 
     495#if defined key_foam_medusa 
     496         WRITE(numout,*)     & 
     497         &   ' key_foam_medusa                                                        = ACTIVE' 
     498#else 
     499         WRITE(numout,*)     & 
     500         &   ' key_foam_medusa                                                        = INACTIVE' 
     501#endif 
    492502         WRITE(numout,*) ' ' 
    493503 
     
    10301040      xthetarem = 0. 
    10311041      xo2min    = 0. 
     1042# if defined key_foam_medusa 
     1043      xobs_xco2a = 0. 
     1044# endif 
    10321045 
    10331046      !READ(numnatm,natroam) 
     
    10491062!!       xthetarem :  oxygen consumption by carbon remineralisation 
    10501063!!       xo2min    :  oxygen minimum concentration 
     1064# if defined key_foam_medusa 
     1065!!       xobs_xco2a : observed atmospheric xCO2 
     1066# endif 
    10511067 
    10521068      IF(lwp) THEN 
     
    10661082          WRITE(numout,*)     & 
    10671083          &   ' oxygen minimum concentration                               xo2min      = ', xo2min 
     1084# if defined key_foam_medusa 
     1085          WRITE(numout,*)     & 
     1086          &   ' observed atmospheric xCO2                                  xobs_xco2a  = ', xobs_xco2a 
     1087# endif 
    10681088       ENDIF 
    10691089 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/trcrst.F90

    r9262 r9292  
    4545   USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable 
    4646   USE trcstat 
     47#if defined key_foam_medusa 
     48   USE obs_const, ONLY: obfillflt  ! Observation operator fill value 
     49#endif 
    4750 
    4851   IMPLICIT NONE 
     
    338341         IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...'  
    339342      ENDIF 
     343      ! 
     344#  if defined key_foam_medusa 
     345      !! 2D fields of pCO2 and fCO2 for observation operator on first timestep 
     346      IF( iom_varid( numrtr, 'PCO2W', ldstop = .FALSE. ) > 0 ) THEN 
     347         IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 present - reading in ...' 
     348         CALL iom_get( numrtr, jpdom_autoglo, 'PCO2W',  f2_pco2w(:,:)  ) 
     349         CALL iom_get( numrtr, jpdom_autoglo, 'FCO2W',  f2_fco2w(:,:)  ) 
     350      ELSE 
     351         IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 absent - setting to fill ...' 
     352         f2_pco2w(:,:) = obfillflt * tmask(:,:,1) 
     353         f2_fco2w(:,:) = obfillflt * tmask(:,:,1) 
     354      ENDIF 
     355#  endif 
    340356# endif 
    341  
     357# if defined key_foam_medusa 
     358      !! Fields for ocean colour assimilation on first timestep 
     359      IF( iom_varid( numrtr, 'pgrow_avg', ldstop = .FALSE. ) > 0 ) THEN 
     360         IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg present - reading in ...' 
     361         CALL iom_get( numrtr, jpdom_autoglo, 'pgrow_avg',  pgrow_avg(:,:)  ) 
     362         CALL iom_get( numrtr, jpdom_autoglo, 'ploss_avg',  ploss_avg(:,:)  ) 
     363         CALL iom_get( numrtr, jpdom_autoglo, 'phyt_avg',   phyt_avg(:,:)   ) 
     364         CALL iom_get( numrtr, jpdom_autoglo, 'mld_max',    mld_max(:,:)    ) 
     365      ELSE 
     366         IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg absent - setting to zero ...' 
     367         pgrow_avg(:,:) = 0.0 
     368         ploss_avg(:,:) = 0.0 
     369         phyt_avg(:,:)  = 0.0 
     370         mld_max(:,:)   = 0.0 
     371      ENDIF 
     372# endif 
    342373 
    343374#endif 
     
    510541      call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG') 
    511542      !! 
     543#  if defined key_foam_medusa 
     544      !! Fields for observation operator on first timestep 
     545      IF(lwp) WRITE(numout,*) ' MEDUSA OBS fields - writing out ...' 
     546      CALL iom_rstput( kt, nitrst, numrtw, 'PCO2W', f2_pco2w(:,:)  ) 
     547      CALL iom_rstput( kt, nitrst, numrtw, 'FCO2W', f2_fco2w(:,:)  ) 
     548#  endif 
     549# endif 
     550# if defined key_foam_medusa 
     551      !! Fields for assimilation on first timestep 
     552      IF(lwp) WRITE(numout,*) ' MEDUSA ASM fields - writing out ...' 
     553      CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg', pgrow_avg(:,:) ) 
     554      CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg', ploss_avg(:,:) ) 
     555      CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg',  phyt_avg(:,:)  ) 
     556      CALL iom_rstput( kt, nitrst, numrtw, 'mld_max',   mld_max(:,:)   ) 
    512557# endif 
    513558!! 
Note: See TracChangeset for help on using the changeset viewer.