Changeset 9381


Ignore:
Timestamp:
2018-03-07T16:19:58+01:00 (2 years ago)
Author:
dford
Message:

Implement pH balancing.

Location:
branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90

    r9333 r9381  
    6262      & asm_pco2_bal 
    6363   USE par_medusa             ! MEDUSA parameters 
    64    USE mocsy_mainmod, ONLY: & ! MEDUSA carbonate system 
     64   USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system 
     65      & mocsy_interface 
     66   USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion 
    6567      & f2pCO2 
    6668   USE sms_medusa, ONLY:    & ! MEDUSA diagnostic variables 
     
    177179#endif 
    178180 
     181#  include "domzgr_substitute.h90" 
     182 
    179183CONTAINS 
    180184 
     
    215219      ENDIF 
    216220 
    217       IF ( ( ln_balwri ).AND.( .NOT. ln_phytobal ).AND. & 
    218          & ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ).AND.( .NOT. ln_pphinc ) ) THEN 
    219          CALL ctl_warn( ' Balancing increments are only calculated for ocean colour', & 
    220             &           ' with ln_phytobal, or for pCO2, fCO2, or pH.',               & 
    221             &           ' Not the case, so ln_balwri will be set to .false.') 
     221      IF ( ( ln_balwri ).AND.                                         & 
     222         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. & 
     223         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. & 
     224         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. & 
     225         & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. & 
     226         & ( .NOT. ln_pchltotinc  ).AND.( .NOT. ln_pphinc      ).AND. & 
     227         & ( .NOT. ln_spco2inc    ).AND.( .NOT. ln_sfco2inc    ) ) THEN 
     228         CALL ctl_warn( ' Balancing increments not being calculated', & 
     229            &           ' so ln_balwri will be set to .false.') 
    222230         ln_balwri = .FALSE. 
    223231      ENDIF 
     
    664672         CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate   ) 
    665673 
    666          IF ( ln_slchltotinc ) THEN 
     674         IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
     675            & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
     676            & ln_slphynoninc ) THEN 
    667677#if defined key_medusa 
    668             CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', phyto2d_balinc(:,:,:,jpchn) ) 
    669             CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', phyto2d_balinc(:,:,:,jpchd) ) 
     678            CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) ) 
     679            CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) ) 
     680            CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) ) 
     681            CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) ) 
     682            CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) ) 
    670683            IF ( ln_phytobal ) THEN 
    671                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', phyto2d_balinc(:,:,:,jpphn) ) 
    672                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', phyto2d_balinc(:,:,:,jpphd) ) 
    673                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', phyto2d_balinc(:,:,:,jppds) ) 
    674                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', phyto2d_balinc(:,:,:,jpzmi) ) 
    675                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', phyto2d_balinc(:,:,:,jpzme) ) 
    676                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', phyto2d_balinc(:,:,:,jpdin) ) 
    677                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', phyto2d_balinc(:,:,:,jpsil) ) 
    678                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', phyto2d_balinc(:,:,:,jpfer) ) 
    679                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', phyto2d_balinc(:,:,:,jpdet) ) 
    680                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', phyto2d_balinc(:,:,:,jpdtc) ) 
    681                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', phyto2d_balinc(:,:,:,jpdic) ) 
    682                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', phyto2d_balinc(:,:,:,jpalk) ) 
    683                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', phyto2d_balinc(:,:,:,jpoxy) ) 
     684               CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) ) 
     685               CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) ) 
     686               CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) ) 
     687               CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) ) 
     688               CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) ) 
     689               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) ) 
     690               CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) ) 
     691               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) ) 
     692               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) ) 
     693               CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) ) 
    684694            ENDIF 
    685695#elif defined key_hadocc 
    686             CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', phyto2d_balinc(:,:,:,jp_had_phy) ) 
     696            CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) ) 
    687697            IF ( ln_phytobal ) THEN 
    688                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', phyto2d_balinc(:,:,:,jp_had_nut) ) 
    689                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) ) 
    690                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', phyto2d_balinc(:,:,:,jp_had_pdn) ) 
    691                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', phyto2d_balinc(:,:,:,jp_had_dic) ) 
    692                CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', phyto2d_balinc(:,:,:,jp_had_alk) ) 
     698               CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) ) 
     699               CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) ) 
     700               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) ) 
     701               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) ) 
     702               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) ) 
    693703            ENDIF 
     704#endif 
     705         ENDIF 
     706 
     707         IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN 
     708#if defined key_medusa 
     709            CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) ) 
     710            CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) ) 
     711            CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) ) 
     712            CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) ) 
     713            CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) ) 
     714#elif defined key_hadocc 
     715            CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) ) 
    694716#endif 
    695717         ENDIF 
     
    697719         IF ( ln_spco2inc ) THEN 
    698720#if defined key_medusa 
    699             CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jpdic) ) 
    700             CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jpalk) ) 
     721            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) ) 
     722            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) ) 
    701723#elif defined key_hadocc 
    702             CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) ) 
    703             CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) ) 
     724            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) ) 
     725            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) ) 
    704726#endif 
    705727         ELSE IF ( ln_sfco2inc ) THEN 
    706728#if defined key_medusa 
    707             CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jpdic) ) 
    708             CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jpalk) ) 
     729            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) ) 
     730            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) ) 
    709731#elif defined key_hadocc 
    710             CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) ) 
    711             CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) ) 
     732            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) ) 
     733            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) ) 
     734#endif 
     735         ENDIF 
     736 
     737         IF ( ln_pphinc ) THEN 
     738#if defined key_medusa 
     739            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) ) 
     740            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) ) 
     741#elif defined key_hadocc 
     742            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) ) 
     743            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) ) 
    712744#endif 
    713745         ENDIF 
     
    944976      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights 
    945977      ! 
    946       INTEGER                          :: ji, jj, jk ! Loop counters 
    947       INTEGER                          :: it         ! Index 
    948       REAL(wp)                         :: zincwgt    ! IAU weight for time step 
    949       REAL(wp)                         :: zfrac_chn  ! Fraction of jpchn 
    950       REAL(wp)                         :: zfrac_chd  ! Fraction of jpchd 
    951       REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc    ! Chlorophyll increments 
    952       REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl    ! Chlorophyll background 
     978      INTEGER                          :: ji, jj, jk   ! Loop counters 
     979      INTEGER                          :: it           ! Index 
     980      REAL(wp)                         :: zincwgt      ! IAU weight for timestep 
     981      REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn 
     982      REAL(wp)                         :: zfrac_chd    ! Fraction of jpchd 
     983      REAL(wp)                         :: zrat_phn_chn ! jpphn:jpchn ratio 
     984      REAL(wp)                         :: zrat_phd_chd ! jpphd:jpchd ratio 
     985      REAL(wp)                         :: zrat_pds_chd ! jppds:jpchd ratio 
     986      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc      ! Chlorophyll increments 
     987      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl      ! Chlorophyll background 
    953988      !!------------------------------------------------------------------------ 
    954989 
     
    10061041                     phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn 
    10071042                     phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd 
     1043                     zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn) 
     1044                     zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd) 
     1045                     zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd) 
     1046                     phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn 
     1047                     phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd 
     1048                     phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd 
    10081049                  ENDIF 
    10091050               END DO 
     
    11921233 
    11931234#if defined key_medusa 
    1194          ! Account for logchl balancing if required 
    1195          IF ( ln_slchltotinc .AND. ln_phytobal ) THEN 
     1235         ! Account for phytoplankton balancing if required 
     1236         IF ( ln_phytobal ) THEN 
    11961237            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic) 
    11971238            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk) 
     
    12061247 
    12071248#elif defined key_hadocc 
    1208          ! Account for slchltot balancing if required 
    1209          IF ( ln_slchltotinc .AND. ln_phytobal ) THEN 
     1249         ! Account for phytoplankton balancing if required 
     1250         IF ( ln_phytobal ) THEN 
    12101251            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic) 
    12111252            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk) 
     
    13791420      !! ** Purpose : Apply the pH assimilation increments. 
    13801421      !! 
    1381       !! ** Method  : Calculate increments to state variables using carbon 
    1382       !!              balancing. 
     1422      !! ** Method  : Calculate increments to DIC and alkalinity from pH 
     1423      !!              Use a similar approach to the pCO2 scheme 
    13831424      !!              Direct initialization or Incremental Analysis Updating. 
    13841425      !! 
    13851426      !! ** Action  :  
    13861427      !!------------------------------------------------------------------------ 
    1387       INTEGER, INTENT(IN)                   :: kt           ! Current time step 
    1388       LOGICAL, INTENT(IN)                   :: ll_asmdin    ! Flag for direct initialisation 
    1389       LOGICAL, INTENT(IN)                   :: ll_asmiau    ! Flag for incremental analysis update 
    1390       INTEGER, INTENT(IN)                   :: kcycper      ! Dimension of pwgtiau 
    1391       REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights 
    1392       LOGICAL, INTENT(IN)                   :: ll_trainc    ! Flag for T/S increments 
     1428      INTEGER,                          INTENT(IN) :: kt        ! Current time step 
     1429      LOGICAL,                          INTENT(IN) :: ll_asmdin ! Flag for direct initialisation 
     1430      LOGICAL,                          INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update 
     1431      INTEGER,                          INTENT(IN) :: kcycper   ! Dimension of pwgtiau 
     1432      REAL(wp), DIMENSION(kcycper),     INTENT(IN) :: pwgtiau   ! IAU weights 
     1433      LOGICAL,                          INTENT(IN) :: ll_trainc ! Flag for T/S increments 
    13931434      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments 
    13941435      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments 
    1395       !!------------------------------------------------------------------------ 
    1396        
    1397       CALL ctl_stop( ' pH balancing not yet implemented' ) 
    1398        
    1399       ! See solve_at_general line 281 of mocsy_phsolvers.F90 
    1400       ! 
    1401       ! Or, call to mocsy_interface line 158 of carb_chem.F90 
    1402       ! Input variables (rest are output) are: 
    1403       ! ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), & 
    1404       ! zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), & 
    1405       ! f_pp0(ji,jj),fsdept(ji,jj,jk),       & 
    1406       ! gphit(ji,jj),f_kw660(ji,jj),         & 
    1407       ! f_xco2a(ji,jj),1 
    1408       ! 
    1409       ! ztmp    = tsn(:,:,:,jp_tem) 
    1410       ! zsal    = tsn(:,:,:,jp_sal) 
    1411       ! zalk    = trn(:,:,:,jpalk) 
    1412       ! zdic    = trn(:,:,:,jpdic) 
    1413       ! zsil    = trn(:,:,:,jpsil) 
    1414       ! zpho    = trn(:,:,:,jpdin) / 16.0 
    1415       ! f_pp0   = 1.0 
    1416       ! fsdept  = fsdept(:,:,:) 
    1417       ! gphit   = gphit(:,:) 
    1418       ! f_kw660 = 1.0 
    1419       ! f_xco2a = f_xco2a(:,:) 
    1420        
     1436       
     1437      REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation 
     1438      REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation 
     1439      REAL(wp)                         :: PH_OUT         ! pH from pH calculation 
     1440      REAL(wp)                         :: ph_bkg         ! pH from pH calculation 
     1441      REAL(wp)                         :: ph_dic         ! pH from pH calculation 
     1442      REAL(wp)                         :: ph_alk         ! pH from pH calculation 
     1443      REAL(wp)                         :: dph_ddic       ! Change in pH wrt DIC 
     1444      REAL(wp)                         :: dph_dalk       ! Change in pH wrt alk 
     1445      REAL(wp)                         :: weight         ! Increment weighting 
     1446      REAL(wp)                         :: zincwgt        ! IAU weight for current time step 
     1447      REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp   ! DIC background 
     1448      REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp   ! Alkalinity background 
     1449      REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp   ! DIN background 
     1450      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp   ! Silicate background 
     1451      REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp   ! Temperature background 
     1452      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp   ! Salinity background 
     1453      REAL(wp), DIMENSION(19)          :: dummy_out      ! Dummy variables 
     1454      INTEGER                          :: ji, jj, jk, jx ! Loop counters 
     1455      INTEGER                          :: it             ! Index 
     1456      !!------------------------------------------------------------------------ 
     1457 
     1458#if ! (defined key_medusa && defined key_foam_medusa) 
     1459      CALL ctl_stop( ' pH balancing only implemented for MEDUSA' ) 
     1460#else 
     1461 
     1462      IF ( kt <= nit000 ) THEN 
     1463 
     1464         !---------------------------------------------------------------------- 
     1465         ! DIC and alkalinity balancing 
     1466         !---------------------------------------------------------------------- 
     1467 
     1468         ! Account for physics assimilation if required 
     1469         IF ( ll_trainc ) THEN 
     1470            tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:) 
     1471            sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:) 
     1472         ELSE 
     1473            tem_bkg_temp(:,:,:) = tsn(:,:,:,1) 
     1474            sal_bkg_temp(:,:,:) = tsn(:,:,:,2) 
     1475         ENDIF 
     1476 
     1477         ! Account for phytoplankton balancing if required 
     1478         IF ( ln_phytobal ) THEN 
     1479            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic) 
     1480            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk) 
     1481            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin) 
     1482            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil) 
     1483         ELSE 
     1484            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) 
     1485            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) 
     1486            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) 
     1487            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) 
     1488         ENDIF 
     1489          
     1490         ! Account for pCO2 balancing if required 
     1491         IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN 
     1492            dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic) 
     1493            alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk) 
     1494         ENDIF 
     1495          
     1496         ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk 
     1497         ! This requires three calls to the MOCSY carbonate package 
     1498         ! One to find pH at background DIC and alkalinity 
     1499         ! One to find pH when DIC is increased by 10 
     1500         ! One to find pH when alkalinity is increased by 10 
     1501         ! Then calculate the gradients 
     1502 
     1503         DO jk = 1, jpk 
     1504            DO jj = 1, jpj 
     1505               DO ji = 1, jpi 
     1506 
     1507                  IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN 
     1508 
     1509                     DO jx = 1, 3 
     1510                        IF ( jx == 1 ) THEN 
     1511                           DIC_IN = dic_bkg_temp(ji,jj,jk) 
     1512                           ALK_IN = alk_bkg_temp(ji,jj,jk) 
     1513                        ELSE IF ( jx == 2 ) THEN 
     1514                           DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch 
     1515                           ALK_IN = alk_bkg_temp(ji,jj,jk) 
     1516                        ELSE IF ( jx == 3 ) THEN 
     1517                           DIC_IN = dic_bkg_temp(ji,jj,jk) 
     1518                           ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch 
     1519                        ENDIF 
     1520 
     1521                        CALL mocsy_interface(tsn(ji,jj,jk,jp_tem),         & ! IN:  temperature 
     1522                           &                 tsn(ji,jj,jk,jp_sal),         & !      salinity 
     1523                           &                 ALK_IN,                       & !      alkalinity 
     1524                           &                 DIC_IN,                       & !      DIC 
     1525                           &                 sil_bkg_temp(ji,jj,jk),       & !      silicate 
     1526                           &                 din_bkg_temp(ji,jj,jk)/16.0,  & !      phosphate (Redfield from DIN) 
     1527                           &                 1.0,                          & !      atmospheric pressure (dummy) 
     1528                           &                 fsdept(ji,jj,jk),             & !      depth 
     1529                           &                 gphit(ji,jj),                 & !      latitude 
     1530                           &                 1.0,                          & !      gas transfer (dummy) 
     1531                           &                 1.0,                          & !      atmospheric xCO2 (dummy) 
     1532                           &                 1,                            & !      number of points 
     1533                           &                 PH_OUT,                       & ! OUT: pH 
     1534                           &                 dummy_out(1),  dummy_out(2),  & !      stuff we don't care about 
     1535                           &                 dummy_out(3),  dummy_out(4),  & 
     1536                           &                 dummy_out(5),  dummy_out(6),  & 
     1537                           &                 dummy_out(7),  dummy_out(8),  & 
     1538                           &                 dummy_out(9),  dummy_out(10), & 
     1539                           &                 dummy_out(11), dummy_out(12), & 
     1540                           &                 dummy_out(13), dummy_out(14), & 
     1541                           &                 dummy_out(15), dummy_out(16), & 
     1542                           &                 dummy_out(17), dummy_out(18), & 
     1543                           &                 dummy_out(19) ) 
     1544 
     1545                        IF ( jx == 1 ) THEN 
     1546                           ph_bkg = PH_OUT 
     1547                        ELSE IF ( jx == 2 ) THEN 
     1548                           ph_dic = PH_OUT 
     1549                        ELSE IF ( jx == 3 ) THEN 
     1550                           ph_alk = PH_OUT 
     1551                        ENDIF 
     1552                     END DO 
     1553 
     1554                     dph_ddic = (ph_dic - ph_bkg) / zsearch 
     1555                     dph_dalk = (ph_alk - ph_bkg) / zsearch 
     1556                     weight   = pph_bkginc(ji,jj,jk) / ( (dph_ddic*dph_ddic) + (dph_dalk*dph_dalk) ) 
     1557 
     1558                     ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic 
     1559                     ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk 
     1560                      
     1561                  ENDIF 
     1562                   
     1563               END DO 
     1564            END DO 
     1565         END DO 
     1566 
     1567      ENDIF 
     1568       
     1569      IF ( ll_asmiau ) THEN 
     1570 
     1571         !-------------------------------------------------------------------- 
     1572         ! Incremental Analysis Updating 
     1573         !-------------------------------------------------------------------- 
     1574 
     1575         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
     1576 
     1577            it = kt - nit000 + 1 
     1578            zincwgt = pwgtiau(it)      ! IAU weight for the current time step 
     1579            ! note this is not a tendency so should not be divided by rdt 
     1580 
     1581            IF(lwp) THEN 
     1582               WRITE(numout,*)  
     1583               WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', & 
     1584                  &  kt,' with IAU weight = ', pwgtiau(it) 
     1585               WRITE(numout,*) '~~~~~~~~~~~~' 
     1586            ENDIF 
     1587 
     1588            ! Update the biogeochemical variables 
     1589            ! Add directly to trn and trb, rather than to tra, because tra gets 
     1590            ! reset to zero at the start of trc_stp, called after this routine 
     1591            DO jk = 1, jpkm1 
     1592               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 
     1593                  &                          ph_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
     1594               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 
     1595                  &                          ph_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
     1596            END DO 
     1597 
     1598            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
     1599            ! which is called at end of model run 
     1600 
     1601         ENDIF 
     1602 
     1603      ELSEIF ( ll_asmdin ) THEN  
     1604 
     1605         !-------------------------------------------------------------------- 
     1606         ! Direct Initialization 
     1607         !-------------------------------------------------------------------- 
     1608          
     1609         IF ( kt == nitdin_r ) THEN 
     1610 
     1611            neuler = 0                    ! Force Euler forward step 
     1612 
     1613            ! Initialize the now fields with the background + increment 
     1614            ! Background currently is what the model is initialised with 
     1615            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', & 
     1616               &           ' Background state is taken from model rather than background file' ) 
     1617            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 
     1618               &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) 
     1619            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 
     1620  
     1621            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
     1622            ! which is called at end of model run 
     1623         ENDIF 
     1624         ! 
     1625      ENDIF 
     1626#endif       
    14211627      ! 
    14221628   END SUBROUTINE ph_asm_inc 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r9333 r9381  
    12391239       
    12401240      ! Surface pCO2/fCO2 next 
    1241       IF( ln_sfco2inc .OR. ln_spco2inc ) THEN 
     1241      IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN 
    12421242         CALL pco2_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, & 
    12431243            &               ln_trainc, t_bkginc, s_bkginc ) 
     
    12451245       
    12461246      ! Profile pH next 
    1247       IF( ln_pphinc ) THEN 
     1247      IF ( ln_pphinc ) THEN 
    12481248         CALL ph_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, & 
    12491249            &             ln_trainc, t_bkginc, s_bkginc ) 
Note: See TracChangeset for help on using the changeset viewer.