Changeset 9182


Ignore:
Timestamp:
2018-01-04T17:24:35+01:00 (3 years ago)
Author:
jpalmier
Message:

JPALM - UKESM #475 — merge with dev_r5518_GO6_Carb_Fail_from_GO6_8356 rev 9070-9133

Location:
branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_9163/NEMOGCM/NEMO
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_9163/NEMOGCM/NEMO/OPA_SRC/stpctl.F90

    r6487 r9182  
    5858      INTEGER  ::   ii, ij, ik              ! temporary integers 
    5959      REAL(wp) ::   zumax, zsmin, zssh2     ! temporary scalars 
     60      REAL(wp) ::   ztmax, ztmin            ! temporary scalars 
    6061      INTEGER, DIMENSION(3) ::   ilocu      !  
    6162      INTEGER, DIMENSION(2) ::   ilocs      !  
     
    1481499500  FORMAT (' kt=',i6,' min SSS: ',1pg11.4,', i j: ',2i5) 
    149150 
     151! ==================================================================================================== 
     152! ==================================================================================================== 
     153      !                                              !AXY (25/10/17) 
     154      !                                              !* Test max/min limits of temperature 
     155      !                                              !  ---------------------------------- 
     156      ztmax =  -5.e0  ! arbitrary  low maximum value 
     157      ztmin = 100.e0  ! arbitrary high minimum value 
     158      DO jj = 2, jpjm1 
     159         DO ji = 2, jpim1 
     160            IF( tmask(ji,jj,1) == 1) THEN 
     161               ztmax = MAX(ztmax,tsn(ji,jj,1,jp_tem))     ! find local maximum 
     162               ztmin = MIN(ztmin,tsn(ji,jj,1,jp_tem))     ! find local minimum 
     163            ENDIF 
     164         END DO 
     165      END DO 
     166      IF( lk_mpp )   CALL mpp_max( ztmax )                ! max over the global domain 
     167      IF( lk_mpp )   CALL mpp_min( ztmin )                ! min over the global domain 
     168      ! 
     169      IF( ztmax > 40.) THEN  ! we've got a problem 
     170         IF (lk_mpp) THEN 
     171            CALL mpp_maxloc ( tsn(:,:,1,jp_tem),tmask(:,:,1), ztmax, ii,ij ) 
     172         ELSE 
     173            ilocs = MAXLOC( tsn(:,:,1,jp_tem), mask = tmask(:,:,1) == 1.e0 ) 
     174            ii = ilocs(1) + nimpp - 1 
     175            ij = ilocs(2) + njmpp - 1 
     176         ENDIF 
     177         ! 
     178         IF(lwp) THEN 
     179            WRITE(numout,*) 'stp_ctl:tracer anomaly: *****    WARNING     *****' 
     180            WRITE(numout,*) 'stp_ctl:tracer anomaly: sea surface temperature > 40C' 
     181            WRITE(numout,9600) kt, ztmax, ii, ij 
     182            WRITE(numout,*) 'stp_ctl:tracer anomaly: ***** END OF WARNING *****' 
     183         ENDIF 
     184      ENDIF 
     185      ! 
     186      IF( ztmin < -3.) THEN  ! we've got a problem 
     187         IF (lk_mpp) THEN 
     188            CALL mpp_minloc ( tsn(:,:,1,jp_tem),tmask(:,:,1), ztmin, ii,ij ) 
     189         ELSE 
     190            ilocs = MINLOC( tsn(:,:,1,jp_tem), mask = tmask(:,:,1) == 1.e0 ) 
     191            ii = ilocs(1) + nimpp - 1 
     192            ij = ilocs(2) + njmpp - 1 
     193         ENDIF 
     194         ! 
     195         IF(lwp) THEN 
     196            WRITE(numout,*) 'stp_ctl:tracer anomaly: *****    WARNING     *****' 
     197            WRITE(numout,*) 'stp_ctl:tracer anomaly: sea surface temperature < -3C' 
     198            WRITE(numout,9700) kt, ztmin, ii, ij 
     199            WRITE(numout,*) 'stp_ctl:tracer anomaly: ***** END OF WARNING *****' 
     200         ENDIF 
     201      ENDIF 
     2029600  FORMAT ('stp_ctl:tracer anomaly: kt=',i6,' max SST: ',f16.10,', i j: ',2i5) 
     2039700  FORMAT ('stp_ctl:tracer anomaly: kt=',i6,' min SST: ',f16.10,', i j: ',2i5) 
     204! ==================================================================================================== 
     205! ==================================================================================================== 
    150206       
    151207      IF( lk_c1d )  RETURN          ! No log file in case of 1D vertical configuration 
  • branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_9163/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90

    r8521 r9182  
    6262# endif 
    6363                                   zchd, zchn, zdin, zsil 
    64       USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, tmask 
     64      USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, tmask, mig, mjg 
    6565# if defined key_iomput 
    6666      USE iom,               ONLY: lk_iomput 
     
    9191      USE trcoxy_medusa,     ONLY: trc_oxy_medusa 
    9292# endif 
     93      USE lib_mpp,           ONLY: ctl_stop 
     94      USE trcstat  
    9395 
    9496   !!* Substitution 
     
    122124 
    123125# if defined key_roam 
     126      !! init 
     127      f_fco2w(:,:)       = 0.0 
     128      f_fco2atm(:,:)     = 0.0 
     129      f_schmidtco2(:,:)  = 0.0 
     130      f_kwco2(:,:)       = 0.0 
     131      f_co2starair(:,:)  = 0.0 
     132      f_dpco2(:,:)       = 0.0 
     133      f_rhosw(:,:)       = 0.0 
     134      f_K0(:,:)          = 0.0 
     135      !! air pressure (atm); ultimately this will use air  
     136      !! pressure at the base of the UKESM1 atmosphere  
     137      !!                                      
     138      f_pp0(:,:)   = 1.0 
     139 
     140 
    124141      !!----------------------------------------------------------- 
    125142      !! Air-sea gas exchange 
     
    134151         DO ji = 2,jpim1 
    135152            !! OPEN wet point IF..THEN loop 
    136             if (tmask(ji,jj,1) == 1) then 
     153            IF (tmask(ji,jj,1) == 1) then 
    137154               IF (lk_oasis) THEN 
    138155                  !! use 2D atm xCO2 from atm coupling 
    139156                  f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj) 
     157                  !!!  
     158                  !!! Jpalm test on atm xCO2 
     159                  IF ( (f_xco2a(ji,jj) > 1500 ).OR.(f_xco2a(ji,jj) < 100 ) ) THEN 
     160                    IF(lwp) THEN  
     161                       WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj),          & 
     162                                       ' -- ji =', mig(ji),' jj = ', mjg(jj) 
     163                       CALL ctl_stop( 'MEDUSA - Air-Sea :', 'unrealistic atm xCO2 ' ) 
     164                    ENDIF  
     165                 ENDIF  
    140166               ENDIF 
    141167               !! 
     
    162188               'air-sea: carb-chem kt = ', kt 
    163189               CALL flush(numout) 
     190               !! JPALM add carb print: 
     191               call trc_rst_dia_stat(f_xco2a(:,:), 'f_xco2a') 
     192               call trc_rst_dia_stat(wndm(:,:), 'wndm') 
     193               call trc_rst_dia_stat(f_kw660(:,:), 'f_kw660') 
     194               call trc_rst_dia_stat(ztmp(:,:), 'ztmp') 
     195               call trc_rst_dia_stat(zsal(:,:), 'zsal') 
     196               call trc_rst_dia_stat(zalk(:,:), 'zalk') 
     197               call trc_rst_dia_stat(zdic(:,:), 'zdic') 
     198               call trc_rst_dia_stat(zsil(:,:), 'zsil') 
     199               call trc_rst_dia_stat(zpho(:,:), 'zpho') 
    164200#   endif 
    165201      DO jj = 2,jpjm1 
    166202         DO ji = 2,jpim1 
    167203            if (tmask(ji,jj,1) == 1) then 
    168                !! air pressure (atm); ultimately this will use air  
    169                !! pressure at the base of the UKESM1 atmosphere  
    170                !!                                      
    171                f_pp0(ji,jj)   = 1.0 
    172                !! 
    173                !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp(ji,jj) 
    174                !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm    =', wndm(ji,jj) 
    175                !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i    =', fr_i(ji,jj) 
    176204               !! 
    177205#  if defined key_axy_carbchem 
    178206#   if defined key_mocsy 
     207               !! Jpalm -- 12-09-2017 -- add extra check after reccurent 
     208               !!          carbonate failure in the coupled run. 
     209               !!          must be associated to air-sea flux or air xCO2...i 
     210               !!          Check MOCSY inputs 
     211               IF ( (zsal(ji,jj) > 75.0 ).OR.(zsal(ji,jj) < 0.0 ) .OR.        & 
     212                    (ztmp(ji,jj) > 50.0 ).OR.(ztmp(ji,jj) < -20.0 ) .OR.      & 
     213                    (zalk(ji,jj) > 35.0E2 ).OR.(zalk(ji,jj) <= 0.0 ) .OR.     & 
     214                    (zdic(ji,jj) > 35.0E2 ).OR.(zdic(ji,jj) <= 0.0 ) .OR.     & 
     215                    (f_kw660(ji,jj) > 1.0E-2 ).OR.(f_kw660(ji,jj) < 0.0 ) ) THEN 
     216                  IF(lwp) THEN  
     217                      WRITE(numout,*) ' surface T = ',ztmp(ji,jj) 
     218                      WRITE(numout,*) ' surface S = ',zsal(ji,jj) 
     219                      WRITE(numout,*) ' surface ALK = ',zalk(ji,jj) 
     220                      WRITE(numout,*) ' surface DIC = ',zdic(ji,jj) 
     221                      WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj) 
     222                      WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)    
     223                      WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj) 
     224                      WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj) 
     225                      WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj) 
     226                      WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj) 
     227                      WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj) 
     228                      WRITE(numout,*) ' MOCSY input: ji =', mig(ji),' jj = ', mjg(jj),  & 
     229                                       ' kt = ', kt  
     230                      WRITE(numout,*) 'MEDUSA - Air-Sea INPUT: unrealistic surface Carb. Chemistry' 
     231                      CALL ctl_stop( 'MEDUSA - Air-Sea INPUT: ',             & 
     232                                     'unrealistic surface Carb. Chemistry -- INPUTS' ) 
     233                  ENDIF      
     234               ENDIF      
    179235               !! 
    180236               !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 
     
    201257               f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000. 
    202258               f_dcf(ji,jj)  = f_rhosw(ji,jj) 
     259               !! Jpalm -- 12-09-2017 -- add extra check after reccurent 
     260               !!          carbonate failure in the coupled run. 
     261               !!          must be associated to air-sea flux or air xCO2...i 
     262               !!          Check MOCSY inputs 
     263              IF ( (f_pco2w(ji,jj) > 1.E4 ).OR.(f_pco2w(ji,jj) < 0.0 ) .OR.     & 
     264                   (f_fco2w(ji,jj) > 1.E4 ).OR.(f_fco2w(ji,jj) < 0.0 ) .OR.     &    
     265                   (f_fco2atm(ji,jj) > 1.E4 ).OR.(f_fco2atm(ji,jj) < 0.0 ) .OR.     & 
     266                   (f_co2flux(ji,jj) > 1.E-2 ).OR.(f_co2flux(ji,jj) < -1.E-2 ) .OR.     & 
     267                   (f_dpco2(ji,jj) > 1.E4 ).OR.(f_dpco2(ji,jj) < -1.E4 ) ) THEN 
     268                 IF(lwp) THEN  
     269                     WRITE(numout,*) ' surface T = ',ztmp(ji,jj) 
     270                     WRITE(numout,*) ' surface S = ',zsal(ji,jj) 
     271                     WRITE(numout,*) ' surface ALK = ',zalk(ji,jj) 
     272                     WRITE(numout,*) ' surface DIC = ',zdic(ji,jj) 
     273                     WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj) 
     274                     WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)    
     275                     WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj) 
     276                     WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj) 
     277                     WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj) 
     278                     WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj) 
     279                     WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj) 
     280                     WRITE(numout,*) ' MOCSY output: ji =', mig(ji),' jj = ', mjg(jj),  & 
     281                                       ' kt = ', kt      
     282                     WRITE(numout,*) 'MEDUSA - Air-Sea OUTPUT: unrealistic surface Carb. Chemistry' 
     283                     CALL ctl_stop( 'MEDUSA - Air-Sea OUTPUT: ',            & 
     284                                    'unrealistic surface Carb. Chemistry -- OUTPUTS' ) 
     285                 ENDIF      
     286              ENDIF      
    203287            ENDIF 
    204288         ENDDO 
    205289      ENDDO 
    206290 
     291#   if defined key_debug_medusa 
     292               !! JPALM add carb print: 
     293               call trc_rst_dia_stat(f_pco2w(:,:), 'f_pco2w') 
     294               call trc_rst_dia_stat(f_fco2w(:,:), 'f_fco2w') 
     295               call trc_rst_dia_stat(f_fco2atm(:,:), 'f_fco2atm') 
     296               call trc_rst_dia_stat(f_schmidtco2(:,:), 'f_schmidtco2') 
     297               call trc_rst_dia_stat(f_kwco2(:,:), 'f_kwco2') 
     298               call trc_rst_dia_stat(f_co2starair(:,:), 'f_co2starair') 
     299               call trc_rst_dia_stat(f_co2flux(:,:), 'f_co2flux') 
     300               call trc_rst_dia_stat(f_dpco2(:,:), 'f_dpco2') 
     301#   endif 
    207302#   else    
    208303 
  • branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_9163/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90

    r8521 r9182  
    167167      fslownflux(:,:) = 0.0  
    168168      fslowcflux(:,:) = 0.0  
     169      !! 
     170      !! JPALM -- 21-09-2017 -- needed to debug air-sea carb 
     171      f_xco2a(:,:)  = 0.0 
     172      f_pco2w(:,:)  = 0.0 
     173      f_ph(:,:)     = 0.0 
     174      f_kw660(:,:)  = 0.0 
     175      ztmp(:,:)  = 0.0 
     176      zsal(:,:)  = 0.0 
     177      zalk(:,:)  = 0.0 
     178      zdic(:,:)  = 0.0 
     179      zsil(:,:)  = 0.0 
     180      zpho(:,:)  = 0.0 
     181      f_co2flux(:,:)  = 0.0  
     182      f_pco2atm(:,:)  = 0.0 
     183      f_h2co3(:,:)    = 0.0 
     184      f_hco3(:,:)     = 0.0 
     185      f_co3(:,:)      = 0.0 
     186      f_omarg(:,:)    = 0.0 
     187      f_omcal(:,:)    = 0.0 
    169188      !! 
    170189      !! AXY (08/08/17): zero slow detritus fluxes 
  • branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_9163/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r8442 r9182  
    2323   !!  -   !  2016-11  (A. Yool)              Updated diags for CMIP6 
    2424   !!  -   !  2017-05  (A. Yool)              Added extra DMS calculation 
     25   !!  -   !  2017-11  (J. Palm, A. Yool)     Diagnose tracer excursions 
    2526   !!---------------------------------------------------------------------- 
    2627   !! 
     
    8182                                            gdepw_0, gdepw_n,               & 
    8283                                            nday_year, nsec_day, nyear,     & 
    83                                             rdt, tmask 
     84                                            rdt, tmask, mig, mjg, nimpp,    & 
     85                                            njmpp  
    8486      USE in_out_manager,             ONLY: lwp, numout, nn_date0 
    8587# if defined key_iomput 
     
    8789# endif 
    8890      USE lbclnk,                     ONLY: lbc_lnk 
    89       USE lib_mpp,                    ONLY: ctl_stop 
     91      USE lib_mpp   
    9092      USE oce,                        ONLY: tsb, tsn 
    9193      USE par_kind,                   ONLY: wp 
     
    115117       
    116118      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90 
     119      PUBLIC   trc_bio_exceptional_fix ! here  
    117120 
    118121   !!* Substitution 
     
    177180      !! temporary variables 
    178181      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4 
     182      !! 
     183      !! T and S check temporary variable 
     184      REAL(wp) ::    sumtsn, tsnavg 
     185      INTEGER  ::    summask 
     186      CHARACTER(40) :: charout, charout2, charout3, charout4, charout5 
    179187      !! 
    180188      !!------------------------------------------------------------------ 
     
    450458 
    451459# if defined key_roam 
     460         !! extra MEDUSA-2 tracers 
    452461         DO jj = 2,jpjm1 
    453462            DO ji = 2,jpim1 
     
    456465                  zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 
    457466                  !! dissolved inorganic carbon 
    458                   zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) 
     467                  zdic(ji,jj) = trn(ji,jj,jk,jpdic) 
    459468                  !! alkalinity 
    460                   zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) 
     469                  zalk(ji,jj) = trn(ji,jj,jk,jpalk) 
    461470                  !! oxygen 
    462471                  zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) 
     
    470479                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem) 
    471480                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal) 
    472                   !! 
    473              !! AXY (28/02/14): check input fields 
    474                   if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 
    475                      IF(lwp) WRITE(numout,*)                                 & 
    476                         ' trc_bio_medusa: T WARNING 2D, ',                   & 
    477                         tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          & 
    478                         ' at (', ji, ',', jj, ',', jk, ') at time', kt 
    479            IF(lwp) WRITE(numout,*)                                 & 
    480                         ' trc_bio_medusa: T SWITCHING 2D, ',                 & 
    481                         tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 
    482                      !! temperatur 
    483                      ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) 
    484                   endif 
    485                   if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 
    486                      IF(lwp) WRITE(numout,*)                                 & 
    487                         ' trc_bio_medusa: S WARNING 2D, ',                   & 
    488                         tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          & 
    489                         ' at (', ji, ',', jj, ',', jk, ') at time', kt 
    490                   endif 
    491481               ENDIF 
    492482            ENDDO 
    493483         ENDDO 
    494484# else 
     485         !! diagnostic MEDUSA-1 detrital carbon tracer 
    495486         DO jj = 2,jpjm1 
    496487            DO ji = 2,jpim1 
    497                if (tmask(ji,jj,jk) == 1) then 
     488               IF (tmask(ji,jj,jk) == 1) THEN 
    498489                  !! implicit detrital carbon 
    499490                  zdtc(ji,jj) = zdet(ji,jj) * xthetad 
     
    502493         ENDDO 
    503494# endif 
     495 
     496# if defined key_roam 
     497         !! --------------------------------------------- 
     498         !! JPALM -- 14-12-2017 -- Here, before any exeptionnal crazy value is 
     499         !!              removed, we want to tell to the Master Node that this  
     500         !!              Exceptionnal value did exist 
     501         !! 
     502         Call trc_bio_check(kt) 
     503 
     504         !!================================================================ 
     505    !! AXY (03/11/17): check input fields 
     506         !! tracer values that exceed thresholds can cause carbonate system 
     507         !! failures when passed to MOCSY; temporary temperature excursions 
     508         !! in recent UKESM0.8 runs trigger such failures but are too short 
     509         !! to have physical consequences; this section checks for such 
     510         !! values and amends them using neighbouring values 
     511         !!  
     512         !! the check on temperature here is also carried out at the end of 
     513         !! each model time step and anomalies are reported in the master 
     514         !! ocean.output file; the error reporting below is strictly local 
     515         !! to the relevant ocean.output_XXXX file so will not be visible 
     516         !! unless all nodes are reporting output 
     517         !!================================================================ 
     518         !! 
     519         DO jj = 2,jpjm1 
     520            DO ji = 2,jpim1 
     521               if (tmask(ji,jj,jk) == 1) then 
     522                  !! the thresholds for the four tracers are ... 
     523                  IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT. 40.0  ) .OR.   & 
     524                       (zsal(ji,jj) .LE.  0.0) .OR. (zsal(ji,jj) .GT. 50.0  ) .OR.   & 
     525                       (zdic(ji,jj) .LE.  0.0) .OR. (zdic(ji,jj) .GT. 4.0E3 ) .OR.   & 
     526                       (zalk(ji,jj) .LE.  0.0) .OR. (zalk(ji,jj) .GT. 4.0E3 ) ) THEN 
     527                     !! 
     528                     !! all tracer values are reported in the event of any excursion 
     529                     write(charout,*)  ' Tmp = ', ztmp(ji,jj) 
     530                     write(charout2,*) ' Sal = ', zsal(ji,jj) 
     531                     write(charout3,*) ' DIC = ', zdic(ji,jj) 
     532                     write(charout4,*) ' Alk = ', zalk(ji,jj) 
     533                     write(charout5,*) mig(ji), mjg(jj), jk, kt  
     534                     IF (lwp) CALL ctl_warn( 'trc_bio_medusa: carbonate chemistry WARNING:',  & 
     535                        TRIM(charout),TRIM(charout2),TRIM(charout3),TRIM(charout4),           &  
     536                        'at i, j, k, kt:', TRIM(charout5) ) 
     537                     !! 
     538                     !! Detect, report and correct tracer excursions 
     539                     IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT.  40.0) )             & 
     540                        CALL trc_bio_exceptional_fix(                                         & 
     541                        tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_tem), tmask(ji-1:ji+1,jj-1:jj+1,jk),    & 
     542                        'Tmp', -3.0, 40.0, ztmp(ji,jj) ) 
     543                     !! 
     544                     IF ( (zsal(ji,jj) .LE. 0.0) .OR. (zsal(ji,jj) .GT.  50.0)  )             & 
     545                        CALL trc_bio_exceptional_fix(                                         & 
     546                        tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_sal), tmask(ji-1:ji+1,jj-1:jj+1,jk),    & 
     547                        'Sal', 1.0, 50.0, zsal(ji,jj) ) 
     548                     !! 
     549                     IF ( (zdic(ji,jj) .LE. 0.0) .OR. (zdic(ji,jj) .GT. 4.0E3)  )             & 
     550                        CALL trc_bio_exceptional_fix(                                         & 
     551                        trn(ji-1:ji+1,jj-1:jj+1,jk,jpdic), tmask(ji-1:ji+1,jj-1:jj+1,jk),     & 
     552                        'DIC', 100.0, 4.0E3, zdic(ji,jj) ) 
     553                     !! 
     554                     IF ( (zalk(ji,jj) .LE. 0.0) .OR. (zalk(ji,jj) .GT. 4.0E3)  )             & 
     555                        CALL trc_bio_exceptional_fix(                                         & 
     556                        trn(ji-1:ji+1,jj-1:jj+1,jk,jpalk), tmask(ji-1:ji+1,jj-1:jj+1,jk),     & 
     557                        'Alk', 100.0, 4.0E3, zalk(ji,jj) ) 
     558                  ENDIF 
     559               ENDIF 
     560            ENDDO 
     561         ENDDO 
     562# endif 
     563 
    504564# if defined key_debug_medusa 
    505565         DO jj = 2,jpjm1 
     
    657717   END SUBROUTINE trc_bio_medusa 
    658718 
     719 
     720 
     721   SUBROUTINE trc_bio_exceptional_fix(tiny_var, tiny_mask, var_nm, mini, maxi, varout) 
     722      !! JPALM (27/10/17): This function is called only when abnormal values that  
     723      !! could break the model's carbonate system are fed to MEDUSA 
     724      !!  
     725      !! The function calculates an average tracer value based on the 3x3 cell 
     726      !! neighbourhood around the abnormal cell, and reports this back 
     727      !!  
     728      !! Tracer array values are not modified, but MEDUSA uses "corrected" values 
     729      !! in its calculations 
     730      !! 
     731      !! temporary variables 
     732      REAL(wp), INTENT( in ), DIMENSION(3,3) :: tiny_var, tiny_mask 
     733      CHARACTER(25), INTENT( in )            :: var_nm 
     734      REAL(wp), INTENT( in )                 :: mini, maxi 
     735      REAL(wp), INTENT( out )                :: varout 
     736      REAL(wp)      :: sumtsn, tsnavg 
     737      INTEGER       :: summask 
     738      CHARACTER(25) :: charout1, charout2 
     739      CHARACTER(60) :: charout3, charout4 
     740      INTEGER       :: ii,ij 
     741     
     742      !! point to the center of the 3*3 zoom-grid, to check around 
     743      ii = 2 
     744      ij = 2 
     745      !! Print surounding values to check if isolated Crazy value or  
     746      !! General error 
     747      IF(lwp) WRITE(numout,*)                                 & 
     748         '----------------------------------------------------------------------' 
     749      IF(lwp) WRITE(numout,*)                                 & 
     750         'trc_bio_medusa: 3x3 neighbourhood surrounding abnormal ', TRIM(var_nm) 
     751      IF(lwp) WRITE(numout,9100)                              & 
     752         3, tiny_var(ii-1,ij+1), tiny_var(ii  ,ij+1), tiny_var(ii+1,ij+1) 
     753      IF(lwp) WRITE(numout,9100)                              & 
     754         2, tiny_var(ii-1,ij  ), tiny_var(ii  ,ij  ), tiny_var(ii+1,ij  ) 
     755      IF(lwp) WRITE(numout,9100)                              & 
     756         1, tiny_var(ii-1,ij-1), tiny_var(ii  ,ij-1), tiny_var(ii+1,ij-1) 
     757      IF(lwp) WRITE(numout,*)                                 & 
     758         'trc_bio_medusa: 3x3 land-sea neighbourhood, tmask' 
     759      IF(lwp) WRITE(numout,9100)                              & 
     760         3, tiny_mask(ii-1,ij+1), tiny_mask(ii  ,ij+1), tiny_mask(ii+1,ij+1) 
     761      IF(lwp) WRITE(numout,9100)                              & 
     762         2, tiny_mask(ii-1,ij  ), tiny_mask(ii  ,ij  ), tiny_mask(ii+1,ij  ) 
     763      IF(lwp) WRITE(numout,9100)                              & 
     764         1, tiny_mask(ii-1,ij-1), tiny_mask(ii  ,ij-1), tiny_mask(ii+1,ij-1) 
     765 
     766      !! Correct out of range values 
     767      sumtsn = ( tiny_mask(ii-1,ij+1) * tiny_var(ii-1,ij+1) ) +  & 
     768               ( tiny_mask(ii  ,ij+1) * tiny_var(ii  ,ij+1) ) +  & 
     769               ( tiny_mask(ii+1,ij+1) * tiny_var(ii+1,ij+1) ) +  & 
     770               ( tiny_mask(ii-1,ij  ) * tiny_var(ii-1,ij  ) ) +  & 
     771               ( tiny_mask(ii+1,ij  ) * tiny_var(ii+1,ij  ) ) +  & 
     772               ( tiny_mask(ii-1,ij-1) * tiny_var(ii-1,ij-1) ) +  & 
     773               ( tiny_mask(ii  ,ij-1) * tiny_var(ii  ,ij-1) ) +  & 
     774               ( tiny_mask(ii+1,ij-1) * tiny_var(ii+1,ij-1) ) 
     775      !! 
     776      summask = tiny_mask(ii-1,ij+1) + tiny_mask(ii  ,ij+1)   +  & 
     777                tiny_mask(ii+1,ij+1) + tiny_mask(ii-1,ij  )   +  & 
     778                tiny_mask(ii+1,ij  ) + tiny_mask(ii-1,ij-1)   +  & 
     779                tiny_mask(ii  ,ij-1) + tiny_mask(ii+1,ij-1) 
     780      !! 
     781      IF ( summask .GT. 0 ) THEN 
     782         tsnavg = ( sumtsn / summask ) 
     783         varout = MAX( MIN( maxi, tsnavg), mini ) 
     784      ELSE    
     785         IF (ztmp(ii,ij) .LT. mini )  varout = mini 
     786         IF (ztmp(ii,ij) .GT. maxi )  varout = maxi 
     787      ENDIF 
     788      !! 
     789      write(charout1,9200) tiny_var(ii,ij) 
     790      write(charout2,9200) varout 
     791      write(charout3,*) ' ', charout1, ' -> ', charout2 
     792      write(charout4,*) ' Tracer: ', trim(var_nm) 
     793      IF(lwp) WRITE(numout,*) 'trc_bio_medusa: ** EXCEPTIONAL VALUE SWITCHING **' 
     794      IF(lwp) WRITE(numout,*) charout4  
     795      IF(lwp) WRITE(numout,*) charout3 
     796      IF(lwp) WRITE(numout,*)                                 & 
     797         '----------------------------------------------------------------------' 
     798      IF(lwp) WRITE(numout,*)                                 & 
     799         ' ' 
     8009100  FORMAT('Row:', i1, '  ', e16.6, ' ', e16.6, ' ', e16.6) 
     8019200  FORMAT(e16.6) 
     802 
     803   END SUBROUTINE trc_bio_exceptional_fix  
     804 
     805   SUBROUTINE trc_bio_check(kt) 
     806      !!----------------------------------- 
     807      !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure 
     808      !!                     problem. The model is now able to correct a local 
     809      !!                     crazy value. but does it silently. 
     810      !!                     We need to spread the word to the master node. we 
     811      !!                     don't want the model  to correct values without telling us 
     812      !!                     This module will tell at least when crazy DIC or 
     813      !!                     ALK appears. 
     814      !!------------------------------------- 
     815      REAL(wp)              :: zmax, zmin    ! temporary scalars 
     816      INTEGER               :: ji,jj         ! dummy loop 
     817      INTEGER               :: ii,ij         ! temporary scalars  
     818      INTEGER, DIMENSION(2) :: ilocs         !  
     819      INTEGER, INTENT( in ) :: kt 
     820      !! 
     821      !!========================== 
     822      !! DIC Check 
     823      zmax =  -5.0  ! arbitrary  low maximum value 
     824      zmin =  4.0E4  ! arbitrary high minimum value 
     825      DO jj = 2, jpjm1 
     826         DO ji = 2,jpim1 
     827            IF( tmask(ji,jj,1) == 1) THEN 
     828               zmax = MAX(zmax,zdic(ji,jj))     ! find local maximum 
     829               zmin = MIN(zmin,zdic(ji,jj))     ! find local minimum 
     830            ENDIF 
     831         END DO 
     832      END DO 
     833      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain 
     834      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain 
     835      ! 
     836      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem 
     837         IF (lk_mpp) THEN 
     838            CALL mpp_maxloc ( zdic(:,:),tmask(:,:,1), zmax, ii,ij ) 
     839         ELSE 
     840            ilocs = MAXLOC( zdic(:,:), mask = tmask(:,:,1) == 1. ) 
     841            ii = ilocs(1) + nimpp - 1 
     842            ij = ilocs(2) + njmpp - 1 
     843         ENDIF 
     844         ! 
     845         IF(lwp) THEN 
     846            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
     847            WRITE(numout,*) 'trc_bio:tracer anomaly: sea surface DIC > 4000 ' 
     848            WRITE(numout,9600) kt, zmax, ii, ij 
     849            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
     850         ENDIF 
     851      ENDIF 
     852      ! 
     853      IF( zmin .LE. 0.0) THEN  ! we've got a problem 
     854         IF (lk_mpp) THEN 
     855            CALL mpp_minloc ( zdic(:,:),tmask(:,:,1), zmin, ii,ij ) 
     856         ELSE 
     857            ilocs = MINLOC( zdic(:,:), mask = tmask(:,:,1) == 1. ) 
     858            ii = ilocs(1) + nimpp - 1 
     859            ij = ilocs(2) + njmpp - 1 
     860         ENDIF 
     861         ! 
     862         IF(lwp) THEN 
     863            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
     864            WRITE(numout,*) 'trc_bio:tracer anomaly: sea surface DIC <= 0 ' 
     865            WRITE(numout,9700) kt, zmin, ii, ij 
     866            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
     867         ENDIF 
     868      ENDIF 
     869      !! 
     870      !!========================== 
     871      !! ALKALINITY Check 
     872      zmax =  -5.0  ! arbitrary  low maximum value 
     873      zmin =  4.0E4  ! arbitrary high minimum value 
     874      DO jj = 2, jpjm1 
     875         DO ji = 2,jpim1 
     876            IF( tmask(ji,jj,1) == 1) THEN 
     877               zmax = MAX(zmax,zalk(ji,jj))     ! find local maximum 
     878               zmin = MIN(zmin,zalk(ji,jj))     ! find local minimum 
     879            ENDIF 
     880         END DO 
     881      END DO 
     882      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain 
     883      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain 
     884      ! 
     885      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem 
     886         IF (lk_mpp) THEN 
     887            CALL mpp_maxloc ( zalk(:,:),tmask(:,:,1), zmax, ii,ij ) 
     888         ELSE 
     889            ilocs = MAXLOC( zalk(:,:), mask = tmask(:,:,1) == 1. ) 
     890            ii = ilocs(1) + nimpp - 1 
     891            ij = ilocs(2) + njmpp - 1 
     892         ENDIF 
     893         ! 
     894         IF(lwp) THEN 
     895            WRITE(numout,*) 'trc_bio:tracer anomaly: *****     WARNING    *****' 
     896            WRITE(numout,*) 'trc_bio:tracer anomaly: sea surface Alkalinity > 4000 ' 
     897            WRITE(numout,9800) kt, zmax, ii, ij 
     898            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
     899         ENDIF 
     900      ENDIF 
     901      ! 
     902      IF( zmin .LE. 0.0) THEN  ! we've got a problem 
     903         IF (lk_mpp) THEN 
     904            CALL mpp_minloc ( zalk(:,:),tmask(:,:,1), zmin, ii,ij ) 
     905         ELSE 
     906            ilocs = MINLOC( zalk(:,:), mask = tmask(:,:,1) == 1. ) 
     907            ii = ilocs(1) + nimpp - 1 
     908            ij = ilocs(2) + njmpp - 1 
     909         ENDIF 
     910         ! 
     911         IF(lwp) THEN 
     912            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
     913            WRITE(numout,*) 'trc_bio:tracer anomaly:  sea surface Alkalinity <= 0 ' 
     914            WRITE(numout,9900) kt, zmin, ii, ij 
     915            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
     916         ENDIF 
     917      ENDIF 
     918 
     919 
     9209600  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j: ',2i5) 
     9219700  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j: ',2i5) 
     9229800  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j: ',2i5) 
     9239900  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j: ',2i5) 
     924 
     925   END SUBROUTINE trc_bio_check 
     926 
     927 
    659928#else 
    660929   !!===================================================================== 
Note: See TracChangeset for help on using the changeset viewer.