Changeset 10302 for branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
- Timestamp:
- 2018-11-13T18:21:16+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
r8442 r10302 23 23 !! - ! 2016-11 (A. Yool) Updated diags for CMIP6 24 24 !! - ! 2017-05 (A. Yool) Added extra DMS calculation 25 !! - ! 2017-11 (J. Palm, A. Yool) Diagnose tracer excursions 25 26 !!---------------------------------------------------------------------- 26 27 !! … … 77 78 zfer, zpds, zphd, zphn, zsil, & 78 79 zzme, zzmi 79 USE dom_oce, ONLY: e3t_0, e3t_n, & 80 gdept_0, gdept_n, & 81 gdepw_0, gdepw_n, & 82 nday_year, nsec_day, nyear, & 83 rdt, tmask 80 USE dom_oce, ONLY: e3t_0, gdept_0, gdepw_0, & 81 nday_year, nsec_day, & 82 nyear, nyear_len, ndastp, & 83 nsec_month, & 84 rdt, tmask, mig, mjg, nimpp, & 85 njmpp 86 #if defined key_vvl 87 USE dom_oce, ONLY: e3t_n, gdept_n, gdepw_n 88 #endif 89 84 90 USE in_out_manager, ONLY: lwp, numout, nn_date0 85 # if defined key_iomput86 91 USE iom, ONLY: lk_iomput 87 # endif88 92 USE lbclnk, ONLY: lbc_lnk 89 USE lib_mpp, ONLY: ctl_stop 90 USE oce, ONLY: tsb, tsn 93 USE lib_mpp, ONLY: mpp_max, mpp_maxloc, & 94 mpp_min, mpp_minloc, & 95 ctl_stop, ctl_warn, lk_mpp 96 USE oce, ONLY: tsb, tsn, PCO2a_in_cpl 91 97 USE par_kind, ONLY: wp 92 98 USE par_medusa, ONLY: jpalk, jpchd, jpchn, jpdet, & … … 98 104 !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm 99 105 USE sbc_oce, ONLY: lk_oasis 100 USE sms_medusa, ONLY: hist_pco2 106 USE sms_medusa, ONLY: hist_pco2, co2_yinit, co2_yend, & 107 # if defined key_roam 108 xobs_xco2a, & 109 # endif 110 pgrow_avg, & 111 ploss_avg, phyt_avg, mld_max, & 112 lk_pi_co2, ln_foam_medusa 101 113 USE trc, ONLY: ln_rsttr, nittrc000, trn 102 114 USE bio_medusa_init_mod, ONLY: bio_medusa_init … … 110 122 USE bio_medusa_diag_slice_mod, ONLY: bio_medusa_diag_slice 111 123 USE bio_medusa_fin_mod, ONLY: bio_medusa_fin 124 USE trcstat, ONLY: trc_rst_dia_stat 112 125 113 126 IMPLICIT NONE … … 115 128 116 129 PUBLIC trc_bio_medusa ! called in trcsms_medusa.F90 130 PUBLIC trc_bio_exceptionnal_fix ! here 117 131 118 132 !!* Substitution … … 176 190 !! 177 191 !! temporary variables 178 REAL(wp) :: fq0,fq1,fq2,fq3,fq4 192 REAL(wp) :: fq3,fq4 193 REAL(wp) :: this_y, this_d, this_s, fyear 194 !! 195 !! T and S check temporary variable 196 REAL(wp) :: sumtsn, tsnavg 197 INTEGER :: summask 198 CHARACTER(40) :: charout, charout2, charout3, charout4, charout5 179 199 !! 180 200 !!------------------------------------------------------------------ … … 283 303 !!------------------------------------------------------------------ 284 304 !! 285 !! what's atmospheric pCO2 doing? (data start in 1859) 286 iyr1 = nyear - 1859 + 1 287 iyr2 = iyr1 + 1 288 if (iyr1 .le. 1) then 289 !! before 1860 290 f_xco2a(:,:) = hist_pco2(1) 291 elseif (iyr2 .ge. 242) then 292 !! after 2099 293 f_xco2a(:,:) = hist_pco2(242) 294 else 295 !! just right 296 fq0 = hist_pco2(iyr1) 297 fq1 = hist_pco2(iyr2) 298 fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0) 299 !! AXY (14/06/12): tweaked to make more sense (and be correct) 300 # if defined key_bs_axy_yrlen 301 !! bugfix: for 360d year with HadGEM2-ES forcing 302 fq3 = (real(nday_year) - 1.0 + fq2) / 360.0 303 # else 304 !! original use of 365 days (not accounting for leap year or 305 !! 360d year) 306 fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 307 # endif 308 fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 309 f_xco2a(:,:) = fq4 310 endif 311 # if defined key_axy_pi_co2 312 !! OCMIP pre-industrial pCO2 313 !! f_xco2a(:,:) = 284.725 !! CMIP5 pre-industrial pCO2 314 f_xco2a = 284.317 !! CMIP6 pre-industrial pCO2 315 # endif 316 !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear =', nyear 317 !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day =', real(nsec_day) 318 !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year) 319 !! AXY (29/01/14): remove surplus diagnostics 320 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0 =', fq0 321 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1 =', fq1 322 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2 =', fq2 323 !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3 =', fq3 324 IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 =', f_xco2a(1,1) 305 IF (lk_oasis) THEN 306 !! xCO2 from coupled 307 IF ( ( kt == nittrc000 ) .AND. lwp ) & 308 WRITE(numout,*) '** MEDUSA Atm xCO2 given by the UM **' 309 f_xco2a(:,:) = PCO2a_in_cpl(:,:) 310 !! Check the xCO2 from the UM is OK 311 !! piece of code moved from air-sea.F90 312 !!--- 313 DO jj = 2,jpjm1 314 DO ji = 2,jpim1 315 !! OPEN wet point IF..THEN loop 316 IF (tmask(ji,jj,1) == 1) then 317 !!! Jpalm test on atm xCO2 318 IF ( (f_xco2a(ji,jj) .GT. 10000.0 ).OR. & 319 (f_xco2a(ji,jj) .LE. 0.0 ) ) THEN 320 IF(lwp) THEN 321 WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj), & 322 ' -- ji =', mig(ji),' jj = ', mjg(jj) 323 ENDIF 324 CALL ctl_stop( 'MEDUSA - trc_bio :', & 325 'unrealistic coupled atm xCO2 ' ) 326 ENDIF 327 ENDIF 328 ENDDO 329 ENDDO 330 !!--- 331 ELSEIF (lk_pi_co2) THEN 332 !! OCMIP pre-industrial xCO2 333 IF ( ( kt == nittrc000 ) .AND. lwp ) & 334 WRITE(numout,*) '** MEDUSA Atm xCO2 fixed to pre-industrial value **' 335 !! f_xco2a(:,:) = 284.725 !! CMIP5 pre-industrial pCO2 336 f_xco2a(:,:) = 284.317 !! CMIP6 pre-industrial pCO2 337 ELSEIF ( xobs_xco2a > 0.0 ) THEN 338 IF(lwp) WRITE(numout,*) ' using observed atm pCO2 = ', xobs_xco2a 339 f_xco2a(:,:) = xobs_xco2a 340 ELSE 341 !! xCO2 from file 342 !! AXY - JPALM new interpolation scheme usinf nyear_len 343 this_y = real(nyear) 344 this_d = real(nday_year) 345 this_s = real(nsec_day) 346 !! 347 fyear = this_y + ((this_d - 1) + (this_s / (60. * 60. * 24.))) / real(nyear_len(1)) 348 !! 349 IF ( ( kt == nittrc000 ) .AND. lwp ) THEN 350 WRITE(numout,*) '** MEDUSA Atm xCO2 from file **' 351 WRITE(numout,*) ' MEDUSA year =', this_y 352 WRITE(numout,*) ' Year length =', real(nyear_len(1)) 353 WRITE(numout,*) ' MEDUSA nday_year =', this_d 354 WRITE(numout,*) ' MEDUSA nsec_day =', this_s 355 ENDIF 356 !! 357 !! different case test 358 IF (fyear .LE. co2_yinit) THEN 359 !! before first record -- pre-industrial value 360 f_xco2a(:,:) = hist_pco2(1) 361 ELSEIF (fyear .GE. co2_yend) THEN 362 !! after last record - continue to use the last value 363 f_xco2a(:,:) = hist_pco2(int(co2_yend - co2_yinit + 1.) ) 364 ELSE 365 !! just right 366 iyr1 = int(fyear - co2_yinit) + 1 367 iyr2 = iyr1 + 1 368 fq3 = fyear - real(iyr1) - co2_yinit + 1. 369 fq4 = ((1 - fq3) * hist_pco2(iyr1)) + (fq3 * hist_pco2(iyr2)) 370 f_xco2a(:,:) = fq4 371 !! 372 IF ( ( kt == nittrc000 ) .AND. lwp ) THEN 373 WRITE(numout,*) ' MEDUSA year1 =', iyr1 374 WRITE(numout,*) ' MEDUSA year2 =', iyr2 375 WRITE(numout,*) ' xCO2 year1 =', hist_pco2(iyr1) 376 WRITE(numout,*) ' xCO2 year2 =', hist_pco2(iyr2) 377 WRITE(numout,*) ' Year2 weight =', fq3 378 ENDIF 379 ENDIF 380 ENDIF 381 382 !! Writing xCO2 in output on start and on the 1st tsp of each month 383 IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR. & 384 ( nsec_month .LE. INT(rdt) ) ) THEN 385 IF ( lwp ) WRITE(numout,*) ' *** Atm xCO2 *** -- kt:', kt, & 386 '; current date:', ndastp 387 call trc_rst_dia_stat(f_xco2a(:,:), 'atm xCO2') 388 ENDIF 325 389 # endif 326 390 … … 348 412 !! x * 30d + 1*rdt(i.e: mod = rdt) 349 413 !! ++ need to pass carb-chem output var through restarts 350 If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR. & 351 ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN 414 !!If ( (kt == nitt8rc000 .AND. .NOT.ln_rsttr) .OR. & 415 !! ( (mod(kt*rdt,2592000.)) == rdt) THEN 416 !!============================= 417 !! (Jpalm -- updated for restartability issues) 418 !! We want this to be start of month or if starting afresh from 419 !! climatology - marc 20/6/17 420 !!If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR. & 421 !! ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN 422 !!============================= 423 !! Jpalm -- 15-02-2018 -- need to change 3D carb-chem call freq again. 424 !! previous call did not work, probably the (86400*mod(nn_date0,100) part 425 !! should not be in... 426 !! now use the NEMO calendar tool : nsec_month to be sure to call 427 !! at the beginning of a new month . 428 !! DAF: For FOAM we want to run daily 429 IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR. & 430 ( nsec_month .LE. INT(rdt) ) .OR. & 431 ( nsec_day .LE. INT(rdt) .AND. ln_foam_medusa ) ) THEN 432 IF ( lwp ) WRITE(numout,*) & 433 ' *** 3D carb chem call *** -- kt:', kt, & 434 '; current date:', ndastp 352 435 !!--------------------------------------------------------------- 353 436 !! Calculate the carbonate chemistry for the whole ocean on the first … … 450 533 451 534 # if defined key_roam 535 !! extra MEDUSA-2 tracers 452 536 DO jj = 2,jpjm1 453 537 DO ji = 2,jpim1 … … 456 540 zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 457 541 !! dissolved inorganic carbon 458 zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))542 zdic(ji,jj) = trn(ji,jj,jk,jpdic) 459 543 !! alkalinity 460 zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))544 zalk(ji,jj) = trn(ji,jj,jk,jpalk) 461 545 !! oxygen 462 546 zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) … … 470 554 ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem) 471 555 zsal(ji,jj) = tsn(ji,jj,jk,jp_sal) 472 !!473 !! AXY (28/02/14): check input fields474 if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then475 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', kt479 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 !! temperatur483 ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)484 endif485 if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then486 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', kt490 endif491 556 ENDIF 492 557 ENDDO 493 558 ENDDO 494 559 # else 560 !! diagnostic MEDUSA-1 detrital carbon tracer 495 561 DO jj = 2,jpjm1 496 562 DO ji = 2,jpim1 497 if (tmask(ji,jj,jk) == 1) then563 IF (tmask(ji,jj,jk) == 1) THEN 498 564 !! implicit detrital carbon 499 565 zdtc(ji,jj) = zdet(ji,jj) * xthetad … … 502 568 ENDDO 503 569 # endif 570 571 # if defined key_roam 572 !! --------------------------------------------- 573 !! JPALM -- 14-12-2017 -- Here, before any exeptionnal crazy value is 574 !! removed, we want to tell to the Master Processor that this 575 !! Exceptionnal value did exist 576 !! 577 Call trc_bio_check(kt, jk) 578 579 !!================================================================ 580 !! AXY (03/11/17): check input fields 581 !! tracer values that exceed thresholds can cause carbonate system 582 !! failures when passed to MOCSY; temporary temperature excursions 583 !! in recent UKESM0.8 runs trigger such failures but are too short 584 !! to have physical consequences; this section checks for such 585 !! values and amends them using neighbouring values 586 !! 587 !! the check on temperature here is also carried out at the end of 588 !! each model time step and anomalies are reported in the master 589 !! ocean.output file; the error reporting below is strictly local 590 !! to the relevant ocean.output_XXXX file so will not be visible 591 !! unless all processors are reporting output 592 !!================================================================ 593 !! 594 DO jj = 2,jpjm1 595 DO ji = 2,jpim1 596 if (tmask(ji,jj,jk) == 1) then 597 !! the thresholds for the four tracers are ... 598 IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT. 40.0 ) .OR. & 599 (zsal(ji,jj) .LE. 0.0) .OR. (zsal(ji,jj) .GT. 50.0 ) .OR. & 600 (zdic(ji,jj) .LE. 0.0) .OR. (zdic(ji,jj) .GT. 4.0E3 ) .OR. & 601 (zalk(ji,jj) .LE. 0.0) .OR. (zalk(ji,jj) .GT. 4.0E3 ) ) THEN 602 !! 603 !! all tracer values are reported in the event of any excursion 604 IF (lwp) THEN 605 WRITE(charout,*) ' Tmp = ', ztmp(ji,jj) 606 WRITE(charout2,*) ' Sal = ', zsal(ji,jj) 607 WRITE(charout3,*) ' DIC = ', zdic(ji,jj) 608 WRITE(charout4,*) ' Alk = ', zalk(ji,jj) 609 WRITE(charout5,*) mig(ji), mjg(jj), jk, kt 610 CALL ctl_warn( 'trc_bio_medusa: carbonate chemistry WARNING:', & 611 TRIM(charout),TRIM(charout2),TRIM(charout3),TRIM(charout4), & 612 'at i, j, k, kt:', TRIM(charout5) ) 613 ENDIF 614 !! 615 !! Detect, report and correct tracer excursions 616 IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT. 40.0) ) & 617 CALL trc_bio_exceptionnal_fix( & 618 tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_tem), tmask(ji-1:ji+1,jj-1:jj+1,jk), & 619 'Tmp', -3.0, 40.0, ztmp(ji,jj) ) 620 !! 621 IF ( (zsal(ji,jj) .LE. 0.0) .OR. (zsal(ji,jj) .GT. 50.0) ) & 622 CALL trc_bio_exceptionnal_fix( & 623 tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_sal), tmask(ji-1:ji+1,jj-1:jj+1,jk), & 624 'Sal', 1.0, 50.0, zsal(ji,jj) ) 625 !! 626 IF ( (zdic(ji,jj) .LE. 0.0) .OR. (zdic(ji,jj) .GT. 4.0E3) ) & 627 CALL trc_bio_exceptionnal_fix( & 628 trn(ji-1:ji+1,jj-1:jj+1,jk,jpdic), tmask(ji-1:ji+1,jj-1:jj+1,jk), & 629 'DIC', 100.0, 4.0E3, zdic(ji,jj) ) 630 !! 631 IF ( (zalk(ji,jj) .LE. 0.0) .OR. (zalk(ji,jj) .GT. 4.0E3) ) & 632 CALL trc_bio_exceptionnal_fix( & 633 trn(ji-1:ji+1,jj-1:jj+1,jk,jpalk), tmask(ji-1:ji+1,jj-1:jj+1,jk), & 634 'Alk', 100.0, 4.0E3, zalk(ji,jj) ) 635 ENDIF 636 ENDIF 637 ENDDO 638 ENDDO 639 # endif 640 504 641 # if defined key_debug_medusa 505 642 DO jj = 2,jpjm1 … … 657 794 END SUBROUTINE trc_bio_medusa 658 795 796 797 798 SUBROUTINE trc_bio_exceptionnal_fix(tiny_var, tiny_mask, var_nm, mini, maxi, varout) 799 !! JPALM (27/10/17): This function is called only when abnormal values that 800 !! could break the model's carbonate system are fed to MEDUSA 801 !! 802 !! The function calculates an average tracer value based on the 3x3 cell 803 !! neighbourhood around the abnormal cell, and reports this back 804 !! 805 !! Tracer array values are not modified, but MEDUSA uses "corrected" values 806 !! in its calculations 807 !! 808 !! temporary variables 809 REAL(wp), INTENT( in ), DIMENSION(3,3) :: tiny_var, tiny_mask 810 CHARACTER(3), INTENT( in ) :: var_nm 811 REAL(wp), INTENT( in ) :: mini, maxi 812 REAL(wp), INTENT( out ) :: varout 813 REAL(wp) :: sumtsn, tsnavg 814 INTEGER :: summask 815 CHARACTER(25) :: charout1, charout2 816 CHARACTER(60) :: charout3, charout4 817 INTEGER :: ii,ij 818 819 !! point to the center of the 3*3 zoom-grid, to check around 820 ii = 2 821 ij = 2 822 !! Print surounding values to check if isolated Crazy value or 823 !! General error 824 IF(lwp) THEN 825 WRITE(numout,*) & 826 '----------------------------------------------------------------------' 827 WRITE(numout,*) & 828 'trc_bio_medusa: 3x3 neighbourhood surrounding abnormal ', TRIM(var_nm) 829 WRITE(numout,9100) & 830 3, tiny_var(ii-1,ij+1), tiny_var(ii ,ij+1), tiny_var(ii+1,ij+1) 831 WRITE(numout,9100) & 832 2, tiny_var(ii-1,ij ), tiny_var(ii ,ij ), tiny_var(ii+1,ij ) 833 WRITE(numout,9100) & 834 1, tiny_var(ii-1,ij-1), tiny_var(ii ,ij-1), tiny_var(ii+1,ij-1) 835 WRITE(numout,*) & 836 'trc_bio_medusa: 3x3 land-sea neighbourhood, tmask' 837 WRITE(numout,9100) & 838 3, tiny_mask(ii-1,ij+1), tiny_mask(ii ,ij+1), tiny_mask(ii+1,ij+1) 839 WRITE(numout,9100) & 840 2, tiny_mask(ii-1,ij ), tiny_mask(ii ,ij ), tiny_mask(ii+1,ij ) 841 WRITE(numout,9100) & 842 1, tiny_mask(ii-1,ij-1), tiny_mask(ii ,ij-1), tiny_mask(ii+1,ij-1) 843 ENDIF 844 !! Correct out of range values 845 sumtsn = ( tiny_mask(ii-1,ij+1) * tiny_var(ii-1,ij+1) ) + & 846 ( tiny_mask(ii ,ij+1) * tiny_var(ii ,ij+1) ) + & 847 ( tiny_mask(ii+1,ij+1) * tiny_var(ii+1,ij+1) ) + & 848 ( tiny_mask(ii-1,ij ) * tiny_var(ii-1,ij ) ) + & 849 ( tiny_mask(ii+1,ij ) * tiny_var(ii+1,ij ) ) + & 850 ( tiny_mask(ii-1,ij-1) * tiny_var(ii-1,ij-1) ) + & 851 ( tiny_mask(ii ,ij-1) * tiny_var(ii ,ij-1) ) + & 852 ( tiny_mask(ii+1,ij-1) * tiny_var(ii+1,ij-1) ) 853 !! 854 summask = tiny_mask(ii-1,ij+1) + tiny_mask(ii ,ij+1) + & 855 tiny_mask(ii+1,ij+1) + tiny_mask(ii-1,ij ) + & 856 tiny_mask(ii+1,ij ) + tiny_mask(ii-1,ij-1) + & 857 tiny_mask(ii ,ij-1) + tiny_mask(ii+1,ij-1) 858 !! 859 IF ( summask .GT. 0 ) THEN 860 tsnavg = ( sumtsn / summask ) 861 varout = MAX( MIN( maxi, tsnavg), mini ) 862 ELSE 863 IF (ztmp(ii,ij) .LT. mini ) varout = mini 864 IF (ztmp(ii,ij) .GT. maxi ) varout = maxi 865 ENDIF 866 !! 867 IF (lwp) THEN 868 WRITE(charout1,9200) tiny_var(ii,ij) 869 WRITE(charout2,9200) varout 870 WRITE(charout3,*) ' ', charout1, ' -> ', charout2 871 WRITE(charout4,*) ' Tracer: ', trim(var_nm) 872 !! 873 WRITE(numout,*) 'trc_bio_medusa: ** EXCEPTIONAL VALUE SWITCHING **' 874 WRITE(numout,*) charout4 875 WRITE(numout,*) charout3 876 WRITE(numout,*) '----------------------------------------------------------------------' 877 WRITE(numout,*) ' ' 878 ENDIF 879 880 9100 FORMAT('Row:', i1, ' ', e16.6, ' ', e16.6, ' ', e16.6) 881 9200 FORMAT(e16.6) 882 883 END SUBROUTINE trc_bio_exceptionnal_fix 884 885 SUBROUTINE trc_bio_check(kt, jk) 886 !!----------------------------------- 887 !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure 888 !! problem. The model is now able to correct a local 889 !! crazy value. but does it silently. 890 !! We need to spread the word to the master processor. we 891 !! don't want the model to correct values without telling us 892 !! This module will tell at least when crazy DIC or 893 !! ALK appears. 894 !!------------------------------------- 895 REAL(wp) :: zmax, zmin ! temporary scalars 896 INTEGER :: ji,jj ! dummy loop 897 INTEGER :: ii,ij ! temporary scalars 898 INTEGER, DIMENSION(2) :: ilocs ! 899 INTEGER, INTENT( in ) :: kt, jk 900 !! 901 !!========================== 902 !! DIC Check 903 zmax = -5.0 ! arbitrary low maximum value 904 zmin = 4.0E4 ! arbitrary high minimum value 905 DO jj = 2, jpjm1 906 DO ji = 2,jpim1 907 IF( tmask(ji,jj,1) == 1) THEN 908 zmax = MAX(zmax,zdic(ji,jj)) ! find local maximum 909 zmin = MIN(zmin,zdic(ji,jj)) ! find local minimum 910 ENDIF 911 END DO 912 END DO 913 IF( lk_mpp ) CALL mpp_max( zmax ) ! max over the global domain 914 IF( lk_mpp ) CALL mpp_min( zmin ) ! min over the global domain 915 ! 916 IF( zmax .GT. 4.0E3) THEN ! we've got a problem 917 IF (lk_mpp) THEN 918 CALL mpp_maxloc ( zdic(:,:),tmask(:,:,1), zmax, ii,ij ) 919 ELSE 920 ilocs = MAXLOC( zdic(:,:), mask = tmask(:,:,1) == 1. ) 921 ii = ilocs(1) + nimpp - 1 922 ij = ilocs(2) + njmpp - 1 923 ENDIF 924 ! 925 IF(lwp) THEN 926 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** WARNING *****' 927 WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration > 4000 ' 928 WRITE(numout,9600) kt, zmax, ii, ij, jk 929 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 930 ENDIF 931 ENDIF 932 ! 933 IF( zmin .LE. 0.0) THEN ! we've got a problem 934 IF (lk_mpp) THEN 935 CALL mpp_minloc ( zdic(:,:),tmask(:,:,1), zmin, ii,ij ) 936 ELSE 937 ilocs = MINLOC( zdic(:,:), mask = tmask(:,:,1) == 1. ) 938 ii = ilocs(1) + nimpp - 1 939 ij = ilocs(2) + njmpp - 1 940 ENDIF 941 ! 942 IF(lwp) THEN 943 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** WARNING *****' 944 WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration <= 0 ' 945 WRITE(numout,9700) kt, zmin, ii, ij, jk 946 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 947 ENDIF 948 ENDIF 949 !! 950 !!========================== 951 !! ALKALINITY Check 952 zmax = -5.0 ! arbitrary low maximum value 953 zmin = 4.0E4 ! arbitrary high minimum value 954 DO jj = 2, jpjm1 955 DO ji = 2,jpim1 956 IF( tmask(ji,jj,1) == 1) THEN 957 zmax = MAX(zmax,zalk(ji,jj)) ! find local maximum 958 zmin = MIN(zmin,zalk(ji,jj)) ! find local minimum 959 ENDIF 960 END DO 961 END DO 962 IF( lk_mpp ) CALL mpp_max( zmax ) ! max over the global domain 963 IF( lk_mpp ) CALL mpp_min( zmin ) ! min over the global domain 964 ! 965 IF( zmax .GT. 4.0E3) THEN ! we've got a problem 966 IF (lk_mpp) THEN 967 CALL mpp_maxloc ( zalk(:,:),tmask(:,:,1), zmax, ii,ij ) 968 ELSE 969 ilocs = MAXLOC( zalk(:,:), mask = tmask(:,:,1) == 1. ) 970 ii = ilocs(1) + nimpp - 1 971 ij = ilocs(2) + njmpp - 1 972 ENDIF 973 ! 974 IF(lwp) THEN 975 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** WARNING *****' 976 WRITE(numout,*) 'trc_bio:tracer anomaly: ALK concentration > 4000 ' 977 WRITE(numout,9800) kt, zmax, ii, ij, jk 978 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 979 ENDIF 980 ENDIF 981 ! 982 IF( zmin .LE. 0.0) THEN ! we've got a problem 983 IF (lk_mpp) THEN 984 CALL mpp_minloc ( zalk(:,:),tmask(:,:,1), zmin, ii,ij ) 985 ELSE 986 ilocs = MINLOC( zalk(:,:), mask = tmask(:,:,1) == 1. ) 987 ii = ilocs(1) + nimpp - 1 988 ij = ilocs(2) + njmpp - 1 989 ENDIF 990 ! 991 IF(lwp) THEN 992 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** WARNING *****' 993 WRITE(numout,*) 'trc_bio:tracer anomaly: ALK concentration <= 0 ' 994 WRITE(numout,9900) kt, zmin, ii, ij, jk 995 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 996 ENDIF 997 ENDIF 998 999 1000 9600 FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j k: ',3i5) 1001 9700 FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j k: ',3i5) 1002 9800 FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j k: ',3i5) 1003 9900 FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j k: ',3i5) 1004 1005 END SUBROUTINE trc_bio_check 1006 1007 659 1008 #else 660 1009 !!===================================================================== … … 663 1012 CONTAINS 664 1013 SUBROUTINE trc_bio_medusa( kt ) ! Empty routine 1014 IMPLICIT NONE 665 1015 INTEGER, INTENT( in ) :: kt 666 1016 WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt
Note: See TracChangeset
for help on using the changeset viewer.
