Changeset 8456


Ignore:
Timestamp:
2017-08-23T14:20:41+02:00 (3 years ago)
Author:
dford
Message:

Add pCO2/fCO2 assimilation.

Location:
branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM
Files:
1 added
5 edited

Legend:

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

    r8428 r8456  
    12951295    ln_logchlinc = .false. !  Logical switch for applying logchl increments 
    12961296    ln_logchlbal = .false. !  Logical switch for applying logchl multivariate balancing 
     1297    ln_fco2inc = .false.   !  Logical switch for applying fCO2 increments 
     1298    ln_pco2inc = .false.   !  Logical switch for applying pCO2 increments 
    12971299    ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing 
    12981300    nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1] 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90

    r8428 r8456  
    7474 
    7575         IF ( ln_logchlinc ) THEN 
    76 #if defined key_medusa && defined key_foam_medusa 
     76#if defined key_medusa 
    7777            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', logchl_balinc(:,:,:,jpchn) ) 
    7878            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', logchl_balinc(:,:,:,jpchd) ) 
     
    104104         ENDIF 
    105105 
     106         IF ( ln_pco2inc ) THEN 
     107#if defined key_medusa 
     108            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jpdic) ) 
     109            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jpalk) ) 
     110#elif defined key_hadocc 
     111            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) ) 
     112            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) ) 
     113#endif 
     114         ELSE IF ( ln_fco2inc ) THEN 
     115#if defined key_medusa 
     116            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jpdic) ) 
     117            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jpalk) ) 
     118#elif defined key_hadocc 
     119            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) ) 
     120            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) ) 
     121#endif 
     122         ENDIF 
     123 
    106124         CALL iom_close( inum ) 
    107125      ELSE 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r8440 r8456  
    2323   !!   seaice_asm_inc : Apply the seaice increment 
    2424   !!   logchl_asm_inc : Apply the logchl increment 
     25   !!   pco2_asm_inc   : Apply the pCO2/fCO2 increment 
    2526   !!---------------------------------------------------------------------- 
    2627   USE wrk_nemo         ! Memory Allocation 
     
    6465   USE asmlogchlbal_medusa, ONLY: & 
    6566      & asm_logchl_bal_medusa 
     67   USE asmpco2bal, ONLY: & 
     68      & asm_pco2_bal 
    6669   USE par_medusa 
     70   USE mocsy_mainmod 
    6771#elif defined key_hadocc 
    6872   USE asmlogchlbal_hadocc, ONLY: & 
    6973      & asm_logchl_bal_hadocc 
     74   USE asmpco2bal, ONLY: & 
     75      & asm_pco2_bal 
    7076   USE par_hadocc 
    7177#endif 
     
    8187   PUBLIC   seaice_asm_inc !: Apply the seaice increment 
    8288   PUBLIC   logchl_asm_inc !: Apply the logchl increment 
     89   PUBLIC   pco2_asm_inc   !: Apply the pCO2/fCO2 increment 
    8390 
    8491#if defined key_asminc 
     
    97104   LOGICAL, PUBLIC :: ln_logchlinc = .FALSE.   !: No log10(chlorophyll) increment 
    98105   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 
    99108   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check 
    100109   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
     
    123132 
    124133   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl 
     134   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pco2_bkginc   !: Increment to background pCO2 
    125135#if defined key_top 
    126136   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 
    127138#endif 
    128139#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 
     
    192203         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    193204         &                 ln_logchlinc, ln_logchlbal,                     & 
     205         &                 ln_pco2inc, ln_fco2inc,                         & 
    194206         &                 ln_asmdin, ln_asmiau,                           & 
    195207         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
     
    229241         WRITE(numout,*) '      Logical switch for applying logchl increments         ln_logchlinc = ', ln_logchlinc 
    230242         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 
    231245         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau 
    232246         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg 
     
    287301      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 
    288302         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. & 
    289          &        ( ln_logchlinc ) )) & 
     303         &        ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) )) & 
    290304         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
    291          &                ' and ln_logchlinc is set to .true.', & 
     305         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc is set to .true.', & 
    292306         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', & 
    293307         &                ' Inconsistent options') 
     
    298312 
    299313      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 
    300          & .AND.( .NOT. ln_logchlinc ) )  & 
     314         & .AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) )  & 
    301315         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
    302          &                ' and ln_logchlinc are set to .false. :', & 
     316         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc are set to .false. :', & 
    303317         &                ' The assimilation increments are not applied') 
    304318 
     
    321335         &                ' the cycle interval') 
    322336 
    323       IF ( ( ln_balwri ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN 
    324          CALL ctl_warn( ' Balancing increments are only calculated for logchl', & 
    325             &           ' Not assimilating logchl, so ln_balwri will be set to .false.') 
     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.') 
    326340         ln_balwri = .FALSE. 
    327341      ENDIF 
    328342 
    329       IF ( ( ln_logchlbal ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN 
     343      IF ( ( ln_logchlbal ).AND.( .NOT. ln_logchlinc ) ) THEN 
    330344         CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', & 
    331345            &           ' Not assimilating logchl, so ln_logchlbal will be set to .false.') 
    332346         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' ) 
    333351      ENDIF 
    334352 
     
    436454#endif 
    437455      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 
    438464      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 
    439          &  .OR.( ln_logchlinc ) ) THEN 
     465         &  .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN 
    440466 
    441467         !-------------------------------------------------------------------- 
     
    518544            ! to allow for differences in masks 
    519545            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 
    520559         ENDIF 
    521560 
     
    703742            CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) ) 
    704743#endif 
    705          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 
    706774 
    707775         CALL iom_close( inum ) 
     
    14541522      ! 
    14551523   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 
    14561767    
    14571768   !!====================================================================== 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r8428 r8456  
    162162                IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH 
    163163                IF( ln_logchlinc ) CALL logchl_asm_inc( nit000 - 1 ) 
     164                IF( ln_fco2inc .OR. ln_pco2inc ) CALL pco2_asm_inc( nit000 - 1 ) 
    164165             ENDIF 
    165166          ENDIF 
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8428 r8456  
    255255      IF( lk_asminc .AND. ln_asmiau .AND. ln_logchlinc ) & 
    256256         &               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 
    257259                         CALL trc_stp( kstp )         ! time-stepping 
    258260#endif 
Note: See TracChangeset for help on using the changeset viewer.