New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9431 for branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 – NEMO

Ignore:
Timestamp:
2018-03-26T16:49:51+02:00 (6 years ago)
Author:
dford
Message:

Allow assimilation of non-log chlorophyll, and prepare for assimilation of PFTs.

File:
1 edited

Legend:

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

    r9382 r9431  
    6161   USE asmpco2bal, ONLY:    & ! pCO2 balancing for MEDUSA 
    6262      & asm_pco2_bal 
     63   USE sms_medusa             ! MEDUSA parameters 
    6364   USE par_medusa             ! MEDUSA parameters 
    6465   USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system 
     
    819820   !!=========================================================================== 
    820821 
     822   SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog ) 
     823      !!------------------------------------------------------------------------ 
     824      !!                    ***  ROUTINE asm_bgc_init_incs  *** 
     825      !! 
     826      !! ** Purpose :   convert log increments to non-log 
     827      !! 
     828      !! ** Method  :   need to account for model background, 
     829      !!                cannot simply do 10^log_inc. Need to: 
     830      !!                1) Add log_inc to log10(background) to get log10(analysis) 
     831      !!                2) Take 10^log10(analysis) to get analysis 
     832      !!                3) Subtract background from analysis to get non-log incs 
     833      !! 
     834      !! ** Action  :   return non-log increments 
     835      !! 
     836      !! References :    
     837      !!------------------------------------------------------------------------ 
     838      !! 
     839      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pbkg        ! Background 
     840      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pinc_log    ! Log incs 
     841      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs 
     842      ! 
     843      INTEGER                                     :: ji, jj      ! Loop counters 
     844      !! 
     845      !!------------------------------------------------------------------------ 
     846 
     847      DO jj = 1, jpj 
     848         DO ji = 1, jpi 
     849            IF ( pbkg(ji,jj) > 0.0 ) THEN 
     850               pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + & 
     851                  &                       pinc_log(ji,jj) )      & 
     852                  &                 - pbkg(ji,jj) 
     853            ELSE 
     854               pinc_nonlog(ji,jj) = 0.0 
     855            ENDIF 
     856         END DO 
     857      END DO 
     858 
     859   END SUBROUTINE asm_bgc_unlog_2d 
     860 
     861   !!=========================================================================== 
     862   !!=========================================================================== 
     863   !!=========================================================================== 
     864 
    821865   SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau ) 
    822866      !!------------------------------------------------------------------------ 
     
    837881      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights 
    838882      ! 
    839       INTEGER  :: jk              ! Loop counter 
    840       INTEGER  :: it              ! Index 
    841       REAL(wp) :: zincwgt         ! IAU weight for current time step 
    842       REAL(wp) :: zincper         ! IAU interval in seconds 
    843       !!------------------------------------------------------------------------ 
    844  
    845       IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
    846          & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
    847          & ln_slphynoninc ) THEN 
    848          CALL ctl_stop( ' No PFT assimilation quite yet' ) 
    849       ENDIF 
     883      INTEGER                      :: jk          ! Loop counter 
     884      INTEGER                      :: it          ! Index 
     885      REAL(wp)                     :: zincwgt     ! IAU weight for current time step 
     886      REAL(wp)                     :: zincper     ! IAU interval in seconds 
     887      REAL(wp), DIMENSION(jpi,jpj) :: zmld        ! Mixed layer depth 
     888      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs 
     889      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg 
     890      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs 
     891      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg 
     892#if defined key_medusa && defined key_foam_medusa 
     893      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs 
     894      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg 
     895      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs 
     896      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg 
     897      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs 
     898      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg 
     899      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs 
     900      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg 
     901#endif 
     902      !!------------------------------------------------------------------------ 
    850903       
    851904      IF ( kt <= nit000 ) THEN 
    852905 
     906         ! Un-log any log increments for passing to balancing routines 
     907         ! Total chlorophyll 
     908         IF ( ln_slchltotinc ) THEN 
     909#if defined key_medusa && defined key_foam_medusa 
     910            zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd) 
     911#elif defined key_hadocc 
     912            zbkg_chltot(:,:) = chl_bkg(:,:,1) 
     913#endif 
     914            CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot ) 
     915         ELSE IF ( ln_schltotinc ) THEN 
     916            zinc_chltot(:,:) = schltot_bkginc(:,:) 
     917         ELSE 
     918            zinc_chltot(:,:) = 0.0 
     919         ENDIF 
     920 
     921#if defined key_medusa && defined key_foam_medusa 
     922         ! Diatom chlorophyll 
     923         IF ( ln_slchldiainc ) THEN 
     924            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd) 
     925            CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia ) 
     926         ELSE 
     927            zinc_chldia(:,:) = 0.0 
     928         ENDIF 
     929#endif 
     930 
     931#if defined key_medusa && defined key_foam_medusa 
     932         ! Non-diatom chlorophyll 
     933         IF ( ln_slchlnoninc ) THEN 
     934            zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn) 
     935            CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon ) 
     936         ELSE 
     937            zinc_chlnon(:,:) = 0.0 
     938         ENDIF 
     939#endif 
     940 
     941         ! Total phytoplankton carbon 
     942         IF ( ln_slphytotinc ) THEN 
     943#if defined key_medusa && defined key_foam_medusa 
     944            zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 
     945#elif defined key_hadocc 
     946            zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 
     947#endif 
     948            CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot ) 
     949         ELSE 
     950            zinc_phytot(:,:) = 0.0 
     951         ENDIF 
     952 
     953#if defined key_medusa && defined key_foam_medusa 
     954         ! Diatom phytoplankton carbon 
     955         IF ( ln_slphydiainc ) THEN 
     956            zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd 
     957            CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia ) 
     958         ELSE 
     959            zinc_phydia(:,:) = 0.0 
     960         ENDIF 
     961#endif 
     962 
     963#if defined key_medusa && defined key_foam_medusa 
     964         ! Non-diatom phytoplankton carbon 
     965         IF ( ln_slphynoninc ) THEN 
     966            zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn 
     967            CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon ) 
     968         ELSE 
     969            zinc_phynon(:,:) = 0.0 
     970         ENDIF 
     971#endif 
     972 
     973         ! Select mixed layer 
     974         IF ( ll_asmdin ) THEN 
     975            CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', & 
     976               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' ) 
     977            zmld(:,:) = mld_max_bkg(:,:) 
     978         ELSE 
     979            SELECT CASE( mld_choice_bgc ) 
     980            CASE ( 1 )                   ! Turbocline/mixing depth [W points] 
     981               zmld(:,:) = hmld(:,:) 
     982            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points] 
     983               zmld(:,:) = hmlp(:,:) 
     984            CASE ( 3 )                   ! Kara MLD [Interpolated] 
     985#if defined key_karaml 
     986               IF ( ln_kara ) THEN 
     987                  zmld(:,:) = hmld_kara(:,:) 
     988               ELSE 
     989                  CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', & 
     990                     &           ' but ln_kara=.false.' ) 
     991               ENDIF 
     992#else 
     993               CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', & 
     994                  &           ' but is not defined' ) 
     995#endif 
     996            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points] 
     997               !zmld(:,:) = hmld_tref(:,:) 
     998               CALL ctl_stop( ' hmld_tref mixed layer requested for LogChl assimilation,', & 
     999                  &           ' but is not available in this version' ) 
     1000            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points] 
     1001               zmld(:,:) = hmlpt(:,:) 
     1002            END SELECT 
     1003         ENDIF 
     1004 
    8531005         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt 
    8541006 
    8551007#if defined key_medusa && defined key_foam_medusa 
    856          CALL asm_logchl_bal_medusa( slchltot_bkginc, zincper, mld_choice_bgc, & 
    857             &                        rn_maxchlinc, ln_phytobal, ll_asmdin,  & 
    858             &                        pgrow_avg_bkg, ploss_avg_bkg,           & 
    859             &                        phyt_avg_bkg, mld_max_bkg,              & 
     1008         CALL asm_logchl_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), & 
     1009            &                        zinc_chltot,                         & 
     1010            &                        ln_slchldiainc,                      & 
     1011            &                        zinc_chldia,                         & 
     1012            &                        ln_slchlnoninc,                      & 
     1013            &                        zinc_chlnon,                         & 
     1014            &                        ln_slphytotinc,                      & 
     1015            &                        zinc_phytot,                         & 
     1016            &                        ln_slphydiainc,                      & 
     1017            &                        zinc_phydia,                         & 
     1018            &                        ln_slphynoninc,                      & 
     1019            &                        zinc_phynon,                         & 
     1020            &                        zincper,                             & 
     1021            &                        rn_maxchlinc, ln_phytobal, zmld,     & 
     1022            &                        pgrow_avg_bkg, ploss_avg_bkg,        & 
     1023            &                        phyt_avg_bkg, mld_max_bkg,           & 
    8601024            &                        tracer_bkg, phyto2d_balinc ) 
    8611025#elif defined key_hadocc 
    862          CALL asm_logchl_bal_hadocc( slchltot_bkginc, zincper, mld_choice_bgc, & 
    863             &                        rn_maxchlinc, ln_phytobal, ll_asmdin,  & 
    864             &                        pgrow_avg_bkg, ploss_avg_bkg,           & 
    865             &                        phyt_avg_bkg, mld_max_bkg,              & 
    866             &                        chl_bkg(:,:,1), cchl_p_bkg(:,:,1),      & 
     1026         CALL asm_logchl_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), & 
     1027            &                        zinc_chltot,                         & 
     1028            &                        ln_slphytotinc,                      & 
     1029            &                        zinc_phytot,                         & 
     1030            &                        zincper,                             & 
     1031            &                        rn_maxchlinc, ln_phytobal, zmld,     & 
     1032            &                        pgrow_avg_bkg, ploss_avg_bkg,        & 
     1033            &                        phyt_avg_bkg, mld_max_bkg,           & 
     1034            &                        cchl_p_bkg(:,:,1),                   & 
    8671035            &                        tracer_bkg, phyto2d_balinc ) 
    8681036#else 
    869          CALL ctl_stop( 'Attempting to assimilate slchltot, ', & 
     1037         CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', & 
    8701038            &           'but not defined a biogeochemical model' ) 
    8711039#endif 
     
    8871055            IF(lwp) THEN 
    8881056               WRITE(numout,*)  
    889                WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', & 
     1057               WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', & 
    8901058                  &  kt,' with IAU weight = ', pwgtiau(it) 
    8911059               WRITE(numout,*) '~~~~~~~~~~~~' 
Note: See TracChangeset for help on using the changeset viewer.