Changeset 9385 for branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
- Timestamp:
- 2018-03-08T12:36:19+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
r9257 r9385 81 81 gdept_0, gdept_n, & 82 82 gdepw_0, gdepw_n, & 83 nday_year, nsec_day, nyear, & 83 nday_year, nsec_day, & 84 nyear, nyear_len, ndastp, & 85 nsec_month, & 84 86 rdt, tmask, mig, mjg, nimpp, & 85 87 njmpp … … 92 94 mpp_min, mpp_minloc, & 93 95 ctl_stop, ctl_warn, lk_mpp 94 USE oce, ONLY: tsb, tsn 96 USE oce, ONLY: tsb, tsn, PCO2a_in_cpl 95 97 USE par_kind, ONLY: wp 96 98 USE par_medusa, ONLY: jpalk, jpchd, jpchn, jpdet, & … … 102 104 !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm 103 105 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 105 108 USE trc, ONLY: ln_rsttr, nittrc000, trn 106 109 USE bio_medusa_init_mod, ONLY: bio_medusa_init … … 114 117 USE bio_medusa_diag_slice_mod, ONLY: bio_medusa_diag_slice 115 118 USE bio_medusa_fin_mod, ONLY: bio_medusa_fin 119 USE trcstat, ONLY: trc_rst_dia_stat 116 120 117 121 IMPLICIT NONE … … 181 185 !! 182 186 !! 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 184 189 !! 185 190 !! T and S check temporary variable … … 293 298 !!------------------------------------------------------------------ 294 299 !! 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 335 381 # endif 336 382 … … 358 404 !! x * 30d + 1*rdt(i.e: mod = rdt) 359 405 !! ++ 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 362 425 !!--------------------------------------------------------------- 363 426 !! Calculate the carbonate chemistry for the whole ocean on the first … … 502 565 !! Exceptionnal value did exist 503 566 !! 504 Call trc_bio_check(kt )567 Call trc_bio_check(kt, jk) 505 568 506 569 !!================================================================ … … 810 873 END SUBROUTINE trc_bio_exceptionnal_fix 811 874 812 SUBROUTINE trc_bio_check(kt )875 SUBROUTINE trc_bio_check(kt, jk) 813 876 !!----------------------------------- 814 877 !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure … … 824 887 INTEGER :: ii,ij ! temporary scalars 825 888 INTEGER, DIMENSION(2) :: ilocs ! 826 INTEGER, INTENT( in ) :: kt 889 INTEGER, INTENT( in ) :: kt, jk 827 890 !! 828 891 !!========================== … … 852 915 IF(lwp) THEN 853 916 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 856 919 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 857 920 ENDIF … … 869 932 IF(lwp) THEN 870 933 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 873 936 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 874 937 ENDIF … … 901 964 IF(lwp) THEN 902 965 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 905 968 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 906 969 ENDIF … … 918 981 IF(lwp) THEN 919 982 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 922 985 WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 923 986 ENDIF … … 925 988 926 989 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)990 9600 FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j k: ',3i5) 991 9700 FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j k: ',3i5) 992 9800 FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j k: ',3i5) 993 9900 FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j k: ',3i5) 931 994 932 995 END SUBROUTINE trc_bio_check
Note: See TracChangeset
for help on using the changeset viewer.