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 9385 for branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 – NEMO

Ignore:
Timestamp:
2018-03-08T12:36:19+01:00 (6 years ago)
Author:
marc
Message:

GMED ticket 380. Merge changes from Julien's branches/NERC/dev_r5518_GO6_CO2_cmip.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r9257 r9385  
    8181                                            gdept_0, gdept_n,               & 
    8282                                            gdepw_0, gdepw_n,               & 
    83                                             nday_year, nsec_day, nyear,     & 
     83                                            nday_year, nsec_day,            & 
     84                                            nyear, nyear_len, ndastp,       & 
     85                                            nsec_month,                     & 
    8486                                            rdt, tmask, mig, mjg, nimpp,    & 
    8587                                            njmpp  
     
    9294                                            mpp_min, mpp_minloc,            & 
    9395                                            ctl_stop, ctl_warn, lk_mpp   
    94       USE oce,                        ONLY: tsb, tsn 
     96      USE oce,                        ONLY: tsb, tsn, PCO2a_in_cpl 
    9597      USE par_kind,                   ONLY: wp 
    9698      USE par_medusa,                 ONLY: jpalk, jpchd, jpchn, jpdet,     & 
     
    102104      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm 
    103105      USE sbc_oce,                    ONLY: lk_oasis 
    104       USE sms_medusa,                 ONLY: hist_pco2 
     106      USE sms_medusa,                 ONLY: hist_pco2, co2_yinit, co2_yend, & 
     107                                            lk_pi_co2 
    105108      USE trc,                        ONLY: ln_rsttr, nittrc000, trn 
    106109      USE bio_medusa_init_mod,        ONLY: bio_medusa_init 
     
    114117      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice 
    115118      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin 
     119      USE trcstat,                    ONLY: trc_rst_dia_stat 
    116120 
    117121      IMPLICIT NONE 
     
    181185      !!  
    182186      !! temporary variables 
    183       REAL(wp) ::    fq0,fq1,fq2,fq3,fq4 
     187      REAL(wp) ::    fq3,fq4 
     188      REAL(wp) ::    this_y, this_d, this_s, fyear 
    184189      !! 
    185190      !! T and S check temporary variable 
     
    293298      !!------------------------------------------------------------------ 
    294299      !! 
    295       !! what's atmospheric pCO2 doing? (data start in 1859) 
    296       iyr1  = nyear - 1859 + 1 
    297       iyr2  = iyr1 + 1 
    298       if (iyr1 .le. 1) then 
    299          !! before 1860 
    300          f_xco2a(:,:) = hist_pco2(1) 
    301       elseif (iyr2 .ge. 242) then 
    302          !! after 2099 
    303          f_xco2a(:,:) = hist_pco2(242) 
    304       else 
    305          !! just right 
    306          fq0 = hist_pco2(iyr1) 
    307          fq1 = hist_pco2(iyr2) 
    308          fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0) 
    309          !! AXY (14/06/12): tweaked to make more sense (and be correct) 
    310 #  if defined key_bs_axy_yrlen 
    311          !! bugfix: for 360d year with HadGEM2-ES forcing 
    312          fq3 = (real(nday_year) - 1.0 + fq2) / 360.0   
    313 #  else 
    314          !! original use of 365 days (not accounting for leap year or  
    315          !! 360d year) 
    316          fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 
    317 #  endif 
    318          fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 
    319          f_xco2a(:,:) = fq4 
    320       endif 
    321 #  if defined key_axy_pi_co2 
    322       !! OCMIP pre-industrial pCO2 
    323       !! f_xco2a(:,:) = 284.725  !! CMIP5 pre-industrial pCO2 
    324       f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2  
    325 #  endif 
    326       !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear 
    327       !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day) 
    328       !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year) 
    329       !! AXY (29/01/14): remove surplus diagnostics 
    330       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0 
    331       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1 
    332       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2 
    333       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3 
    334       IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1) 
     300      IF (lk_oasis) THEN  
     301         !! xCO2 from coupled 
     302         IF ( ( kt == nittrc000 ) .AND. lwp )     & 
     303             WRITE(numout,*) '** MEDUSA Atm xCO2 given by the UM **' 
     304         f_xco2a(:,:) = PCO2a_in_cpl(:,:) 
     305         !! Check the xCO2 from the UM is OK 
     306         !! piece of code moved from air-sea.F90 
     307         !!--- 
     308         DO jj = 2,jpjm1 
     309            DO ji = 2,jpim1 
     310               !! OPEN wet point IF..THEN loop 
     311               IF (tmask(ji,jj,1) == 1) then 
     312                  !!! Jpalm test on atm xCO2 
     313                  IF ( (f_xco2a(ji,jj) .GT. 10000.0 ).OR.     & 
     314                       (f_xco2a(ji,jj) .LE. 0.0 ) ) THEN 
     315                     IF(lwp) THEN 
     316                        WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj),       & 
     317                                        ' -- ji =', mig(ji),' jj = ', mjg(jj) 
     318                     ENDIF 
     319                     CALL ctl_stop( 'MEDUSA - trc_bio :',     & 
     320                                    'unrealistic coupled atm xCO2 ' ) 
     321                  ENDIF 
     322               ENDIF 
     323            ENDDO 
     324         ENDDO 
     325         !!---  
     326      ELSEIF (lk_pi_co2) THEN 
     327         !! OCMIP pre-industrial xCO2 
     328         IF ( ( kt == nittrc000 ) .AND. lwp )     & 
     329             WRITE(numout,*) '** MEDUSA Atm xCO2 fixed to pre-industrial value **' 
     330         !! f_xco2a(:,:) = 284.725       !! CMIP5 pre-industrial pCO2 
     331         f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2  
     332      ELSE       
     333         !! xCO2 from file 
     334         !! AXY - JPALM new interpolation scheme usinf nyear_len 
     335         this_y = real(nyear) 
     336         this_d = real(nday_year) 
     337         this_s = real(nsec_day) 
     338         !! 
     339         fyear = this_y + ((this_d - 1) + (this_s / (60. * 60. * 24.))) / real(nyear_len(1)) 
     340         !! 
     341         IF ( ( kt == nittrc000 ) .AND. lwp ) THEN  
     342            WRITE(numout,*) '** MEDUSA Atm xCO2 from file **' 
     343            WRITE(numout,*) ' MEDUSA year      =', this_y 
     344            WRITE(numout,*) ' Year length      =', real(nyear_len(1)) 
     345            WRITE(numout,*) ' MEDUSA nday_year =', this_d 
     346            WRITE(numout,*) ' MEDUSA nsec_day  =', this_s 
     347         ENDIF 
     348         !!  
     349         !! different case test  
     350         IF (fyear .LE. co2_yinit) THEN 
     351            !! before first record -- pre-industrial value 
     352            f_xco2a(:,:) = hist_pco2(1)   
     353         ELSEIF (fyear .GE. co2_yend) THEN 
     354            !! after last record - continue to use the last value 
     355            f_xco2a(:,:) = hist_pco2(int(co2_yend - co2_yinit + 1.) ) 
     356         ELSE 
     357            !! just right 
     358            iyr1 = int(fyear - co2_yinit) + 1 
     359            iyr2 = iyr1 + 1  
     360            fq3 = fyear - real(iyr1) - co2_yinit + 1. 
     361            fq4 = ((1 - fq3) * hist_pco2(iyr1)) + (fq3 * hist_pco2(iyr2)) 
     362            f_xco2a(:,:) = fq4 
     363            !! 
     364            IF ( ( kt == nittrc000 ) .AND. lwp ) THEN  
     365               WRITE(numout,*) ' MEDUSA year1     =', iyr1 
     366               WRITE(numout,*) ' MEDUSA year2     =', iyr2 
     367               WRITE(numout,*) ' xCO2 year1       =', hist_pco2(iyr1) 
     368               WRITE(numout,*) ' xCO2 year2       =', hist_pco2(iyr2) 
     369               WRITE(numout,*) ' Year2 weight     =', fq3 
     370            ENDIF 
     371         ENDIF 
     372      ENDIF 
     373  
     374      !! Writing xCO2 in output on start and on the 1st tsp of each  month 
     375      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
     376           ( nsec_month .LE. INT(rdt) ) )  THEN 
     377           IF ( lwp )  WRITE(numout,*) ' *** Atm xCO2 *** -- kt:', kt,       & 
     378                                       ';  current date:', ndastp 
     379          call trc_rst_dia_stat(f_xco2a(:,:), 'atm xCO2') 
     380      ENDIF 
    335381# endif 
    336382 
     
    358404      !!          x * 30d + 1*rdt(i.e: mod = rdt)    
    359405      !!          ++ need to pass carb-chem output var through restarts 
    360       If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
    361            ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN 
     406      !!If ( (kt == nitt8rc000 .AND. .NOT.ln_rsttr) .OR.          & 
     407      !!     ( (mod(kt*rdt,2592000.)) == rdt) THEN 
     408      !!============================= 
     409      !! (Jpalm -- updated for restartability issues) 
     410      !! We want this to be start of month or if starting afresh from   
     411      !! climatology - marc 20/6/17  
     412      !!If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                         &  
     413      !!     ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN  
     414      !!============================= 
     415      !! Jpalm -- 15-02-2018 -- need to change 3D carb-chem call freq again. 
     416      !! previous call did not work, probably the (86400*mod(nn_date0,100) part 
     417      !! should not be in... 
     418      !! now use the NEMO calendar tool : nsec_month to be sure to call  
     419      !! at the beginning of a new month . 
     420      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
     421           ( nsec_month .LE. INT(rdt) ) )  THEN 
     422           IF ( lwp )  WRITE(numout,*)                                       &          
     423                              ' *** 3D carb chem call *** -- kt:', kt,       & 
     424                              ';  current date:', ndastp  
    362425         !!--------------------------------------------------------------- 
    363426         !! Calculate the carbonate chemistry for the whole ocean on the first 
     
    502565         !!              Exceptionnal value did exist 
    503566         !! 
    504          Call trc_bio_check(kt) 
     567         Call trc_bio_check(kt, jk) 
    505568 
    506569         !!================================================================ 
     
    810873   END SUBROUTINE trc_bio_exceptionnal_fix  
    811874 
    812    SUBROUTINE trc_bio_check(kt) 
     875   SUBROUTINE trc_bio_check(kt, jk) 
    813876      !!----------------------------------- 
    814877      !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure 
     
    824887      INTEGER               :: ii,ij         ! temporary scalars  
    825888      INTEGER, DIMENSION(2) :: ilocs         !  
    826       INTEGER, INTENT( in ) :: kt 
     889      INTEGER, INTENT( in ) :: kt, jk 
    827890      !! 
    828891      !!========================== 
     
    852915         IF(lwp) THEN 
    853916            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
    854             WRITE(numout,*) 'trc_bio:tracer anomaly: sea surface DIC > 4000 ' 
    855             WRITE(numout,9600) kt, zmax, ii, ij 
     917            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration > 4000 ' 
     918            WRITE(numout,9600) kt, zmax, ii, ij, jk 
    856919            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
    857920         ENDIF 
     
    869932         IF(lwp) THEN 
    870933            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
    871             WRITE(numout,*) 'trc_bio:tracer anomaly: sea surface DIC <= 0 ' 
    872             WRITE(numout,9700) kt, zmin, ii, ij 
     934            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration <= 0 ' 
     935            WRITE(numout,9700) kt, zmin, ii, ij, jk 
    873936            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
    874937         ENDIF 
     
    901964         IF(lwp) THEN 
    902965            WRITE(numout,*) 'trc_bio:tracer anomaly: *****     WARNING    *****' 
    903             WRITE(numout,*) 'trc_bio:tracer anomaly: sea surface Alkalinity > 4000 ' 
    904             WRITE(numout,9800) kt, zmax, ii, ij 
     966            WRITE(numout,*) 'trc_bio:tracer anomaly: ALK concentration > 4000 ' 
     967            WRITE(numout,9800) kt, zmax, ii, ij, jk 
    905968            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
    906969         ENDIF 
     
    918981         IF(lwp) THEN 
    919982            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
    920             WRITE(numout,*) 'trc_bio:tracer anomaly:  sea surface Alkalinity <= 0 ' 
    921             WRITE(numout,9900) kt, zmin, ii, ij 
     983            WRITE(numout,*) 'trc_bio:tracer anomaly:  ALK concentration <= 0 ' 
     984            WRITE(numout,9900) kt, zmin, ii, ij, jk 
    922985            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
    923986         ENDIF 
     
    925988 
    926989 
    927 9600  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j: ',2i5) 
    928 9700  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j: ',2i5) 
    929 9800  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j: ',2i5) 
    930 9900  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j: ',2i5) 
     9909600  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j k: ',3i5) 
     9919700  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j k: ',3i5) 
     9929800  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j k: ',3i5) 
     9939900  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j k: ',3i5) 
    931994 
    932995   END SUBROUTINE trc_bio_check 
Note: See TracChangeset for help on using the changeset viewer.