- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/TOP/PISCES/P4Z/p4zche.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/PISCES/P4Z/p4zche.F90
r10425 r13463 130 130 INTEGER :: niter_atgen = jp_maxniter_atgen 131 131 132 !! * Substitutions 133 # include "do_loop_substitute.h90" 134 # include "domzgr_substitute.h90" 132 135 !!---------------------------------------------------------------------- 133 136 !! NEMO/TOP 4.0 , NEMO Consortium (2018) … … 137 140 CONTAINS 138 141 139 SUBROUTINE p4z_che 142 SUBROUTINE p4z_che( Kbb, Kmm ) 140 143 !!--------------------------------------------------------------------- 141 144 !! *** ROUTINE p4z_che *** … … 145 148 !! ** Method : - ... 146 149 !!--------------------------------------------------------------------- 150 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices 147 151 INTEGER :: ji, jj, jk 148 152 REAL(wp) :: ztkel, ztkel1, zt , zsal , zsal2 , zbuf1 , zbuf2 … … 164 168 ! ------------------------------------------------------------- 165 169 IF (neos == -1) THEN 166 salinprac(:,:,:) = ts n(:,:,:,jp_sal) * 35.0 / 35.16504170 salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm) * 35.0 / 35.16504 167 171 ELSE 168 salinprac(:,:,:) = ts n(:,:,:,jp_sal)172 salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm) 169 173 ENDIF 170 174 … … 175 179 ! 0.04°C relative to an exact computation 176 180 ! --------------------------------------------------------------------- 177 DO jk = 1, jpk 178 DO jj = 1, jpj 179 DO ji = 1, jpi 180 zpres = gdept_n(ji,jj,jk) / 1000. 181 za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (salinprac(ji,jj,jk) - 35.0) ) 182 za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 ) 183 tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2 184 END DO 185 END DO 186 END DO 181 DO_3D( 1, 1, 1, 1, 1, jpk ) 182 zpres = gdept(ji,jj,jk,Kmm) / 1000. 183 za1 = 0.04 * ( 1.0 + 0.185 * ts(ji,jj,jk,jp_tem,Kmm) + 0.035 * (salinprac(ji,jj,jk) - 35.0) ) 184 za2 = 0.0075 * ( 1.0 - ts(ji,jj,jk,jp_tem,Kmm) / 30.0 ) 185 tempis(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) - za1 * zpres + za2 * zpres**2 186 END_3D 187 187 ! 188 188 ! CHEMICAL CONSTANTS - SURFACE LAYER … … 245 245 zplat = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) ) 246 246 zc1 = 5.92E-3 + zplat**2 * 5.25E-3 247 zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept _n(ji,jj,jk)))) / 4.42E-6247 zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept(ji,jj,jk,Kmm)))) / 4.42E-6 248 248 zpres = zpres / 10.0 249 249 … … 448 448 END SUBROUTINE p4z_che 449 449 450 SUBROUTINE ahini_for_at(p_hini )450 SUBROUTINE ahini_for_at(p_hini, Kbb ) 451 451 !!--------------------------------------------------------------------- 452 452 !! *** ROUTINE ahini_for_at *** … … 462 462 !!--------------------------------------------------------------------- 463 463 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_hini 464 INTEGER, INTENT(in) :: Kbb ! time level indices 464 465 INTEGER :: ji, jj, jk 465 466 REAL(wp) :: zca1, zba1 … … 471 472 IF( ln_timing ) CALL timing_start('ahini_for_at') 472 473 ! 473 DO jk = 1, jpk 474 DO jj = 1, jpj 475 DO ji = 1, jpi 476 p_alkcb = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 477 p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn) 478 p_bortot = borat(ji,jj,jk) 479 IF (p_alkcb <= 0.) THEN 480 p_hini(ji,jj,jk) = 1.e-3 481 ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN 482 p_hini(ji,jj,jk) = 1.e-10_wp 474 DO_3D( 1, 1, 1, 1, 1, jpk ) 475 p_alkcb = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn) 476 p_dictot = tr(ji,jj,jk,jpdic,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn) 477 p_bortot = borat(ji,jj,jk) 478 IF (p_alkcb <= 0.) THEN 479 p_hini(ji,jj,jk) = 1.e-3 480 ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN 481 p_hini(ji,jj,jk) = 1.e-10_wp 482 ELSE 483 zca1 = p_dictot/( p_alkcb + rtrn ) 484 zba1 = p_bortot/ (p_alkcb + rtrn ) 485 ! Coefficients of the cubic polynomial 486 za2 = aKb3(ji,jj,jk)*(1. - zba1) + ak13(ji,jj,jk)*(1.-zca1) 487 za1 = ak13(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - zca1) & 488 & + ak13(ji,jj,jk)*ak23(ji,jj,jk)*(1. - (zca1+zca1)) 489 za0 = ak13(ji,jj,jk)*ak23(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1)) 490 ! Taylor expansion around the minimum 491 zd = za2*za2 - 3.*za1 ! Discriminant of the quadratic equation 492 ! for the minimum close to the root 493 494 IF(zd > 0.) THEN ! If the discriminant is positive 495 zsqrtd = SQRT(zd) 496 IF(za2 < 0) THEN 497 zhmin = (-za2 + zsqrtd)/3. 483 498 ELSE 484 zca1 = p_dictot/( p_alkcb + rtrn ) 485 zba1 = p_bortot/ (p_alkcb + rtrn ) 486 ! Coefficients of the cubic polynomial 487 za2 = aKb3(ji,jj,jk)*(1. - zba1) + ak13(ji,jj,jk)*(1.-zca1) 488 za1 = ak13(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - zca1) & 489 & + ak13(ji,jj,jk)*ak23(ji,jj,jk)*(1. - (zca1+zca1)) 490 za0 = ak13(ji,jj,jk)*ak23(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1)) 491 ! Taylor expansion around the minimum 492 zd = za2*za2 - 3.*za1 ! Discriminant of the quadratic equation 493 ! for the minimum close to the root 494 495 IF(zd > 0.) THEN ! If the discriminant is positive 496 zsqrtd = SQRT(zd) 497 IF(za2 < 0) THEN 498 zhmin = (-za2 + zsqrtd)/3. 499 ELSE 500 zhmin = -za1/(za2 + zsqrtd) 501 ENDIF 502 p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) 503 ELSE 504 p_hini(ji,jj,jk) = 1.e-7 505 ENDIF 506 ! 507 ENDIF 508 END DO 509 END DO 510 END DO 499 zhmin = -za1/(za2 + zsqrtd) 500 ENDIF 501 p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) 502 ELSE 503 p_hini(ji,jj,jk) = 1.e-7 504 ENDIF 505 ! 506 ENDIF 507 END_3D 511 508 ! 512 509 IF( ln_timing ) CALL timing_stop('ahini_for_at') … … 516 513 !=============================================================================== 517 514 518 SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup )515 SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup, Kbb ) 519 516 520 517 ! Subroutine returns the lower and upper bounds of "non-water-selfionization" … … 525 522 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf 526 523 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup 527 528 p_alknw_inf(:,:,:) = -trb(:,:,:,jppo4) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:) & 524 INTEGER, INTENT(in) :: Kbb ! time level indices 525 526 p_alknw_inf(:,:,:) = -tr(:,:,:,jppo4,Kbb) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:) & 529 527 & - fluorid(:,:,:) 530 p_alknw_sup(:,:,:) = (2. * tr b(:,:,:,jpdic) + 2. * trb(:,:,:,jppo4) + trb(:,:,:,jpsil) ) &528 p_alknw_sup(:,:,:) = (2. * tr(:,:,:,jpdic,Kbb) + 2. * tr(:,:,:,jppo4,Kbb) + tr(:,:,:,jpsil,Kbb) ) & 531 529 & * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:) 532 530 … … 534 532 535 533 536 SUBROUTINE solve_at_general( p_hini, zhi )534 SUBROUTINE solve_at_general( p_hini, zhi, Kbb ) 537 535 538 536 ! Universal pH solver that converges from any given initial value, … … 543 541 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: p_hini 544 542 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: zhi 543 INTEGER, INTENT(in) :: Kbb ! time level indices 545 544 546 545 ! Local variables … … 565 564 IF( ln_timing ) CALL timing_start('solve_at_general') 566 565 567 CALL anw_infsup( zalknw_inf, zalknw_sup )566 CALL anw_infsup( zalknw_inf, zalknw_sup, Kbb ) 568 567 569 568 rmask(:,:,:) = tmask(:,:,:) … … 571 570 572 571 ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree 573 DO jk = 1, jpk 574 DO jj = 1, jpj 575 DO ji = 1, jpi 576 IF (rmask(ji,jj,jk) == 1.) THEN 577 p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 578 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 579 zh_ini = p_hini(ji,jj,jk) 580 581 zdelta = (p_alktot-zalknw_inf(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 582 583 IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN 584 zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) ) 585 ELSE 586 zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2. 587 ENDIF 588 589 zdelta = (p_alktot-zalknw_sup(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 590 591 IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN 592 zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2. 593 ELSE 594 zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) ) 595 ENDIF 596 597 zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk)) 572 DO_3D( 1, 1, 1, 1, 1, jpk ) 573 IF (rmask(ji,jj,jk) == 1.) THEN 574 p_alktot = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn) 575 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 576 zh_ini = p_hini(ji,jj,jk) 577 578 zdelta = (p_alktot-zalknw_inf(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 579 580 IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN 581 zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) ) 582 ELSE 583 zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2. 584 ENDIF 585 586 zdelta = (p_alktot-zalknw_sup(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 587 588 IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN 589 zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2. 590 ELSE 591 zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) ) 592 ENDIF 593 594 zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk)) 595 ENDIF 596 END_3D 597 598 zeqn_absmin(:,:,:) = HUGE(1._wp) 599 600 DO jn = 1, jp_maxniter_atgen 601 DO_3D( 1, 1, 1, 1, 1, jpk ) 602 IF (rmask(ji,jj,jk) == 1.) THEN 603 zfact = rhop(ji,jj,jk) / 1000. + rtrn 604 p_alktot = tr(ji,jj,jk,jptal,Kbb) / zfact 605 zdic = tr(ji,jj,jk,jpdic,Kbb) / zfact 606 zbot = borat(ji,jj,jk) 607 zpt = tr(ji,jj,jk,jppo4,Kbb) / zfact * po4r 608 zsit = tr(ji,jj,jk,jpsil,Kbb) / zfact 609 zst = sulfat (ji,jj,jk) 610 zft = fluorid(ji,jj,jk) 611 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 612 zh = zhi(ji,jj,jk) 613 zh_prev = zh 614 615 ! H2CO3 - HCO3 - CO3 : n=2, m=0 616 znumer_dic = 2.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk) 617 zdenom_dic = ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*(ak13(ji,jj,jk) + zh) 618 zalk_dic = zdic * (znumer_dic/zdenom_dic) 619 zdnumer_dic = ak13(ji,jj,jk)*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh & 620 *(4.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)) 621 zdalk_dic = -zdic*(zdnumer_dic/zdenom_dic**2) 622 623 624 ! B(OH)3 - B(OH)4 : n=1, m=0 625 znumer_bor = akb3(ji,jj,jk) 626 zdenom_bor = akb3(ji,jj,jk) + zh 627 zalk_bor = zbot * (znumer_bor/zdenom_bor) 628 zdnumer_bor = akb3(ji,jj,jk) 629 zdalk_bor = -zbot*(zdnumer_bor/zdenom_bor**2) 630 631 632 ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 633 znumer_po4 = 3.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 634 & + zh*(2.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh* ak1p3(ji,jj,jk)) 635 zdenom_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 636 & + zh*( ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh)) 637 zalk_po4 = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1 638 zdnumer_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 639 & + zh*(4.*ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 640 & + zh*(9.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 641 & + ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) & 642 & + zh*(4.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh * ak1p3(ji,jj,jk) ) ) ) 643 zdalk_po4 = -zpt * (zdnumer_po4/zdenom_po4**2) 644 645 ! H4SiO4 - H3SiO4 : n=1, m=0 646 znumer_sil = aksi3(ji,jj,jk) 647 zdenom_sil = aksi3(ji,jj,jk) + zh 648 zalk_sil = zsit * (znumer_sil/zdenom_sil) 649 zdnumer_sil = aksi3(ji,jj,jk) 650 zdalk_sil = -zsit * (zdnumer_sil/zdenom_sil**2) 651 652 ! HSO4 - SO4 : n=1, m=1 653 aphscale = 1.0 + zst/aks3(ji,jj,jk) 654 znumer_so4 = aks3(ji,jj,jk) * aphscale 655 zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh 656 zalk_so4 = zst * (znumer_so4/zdenom_so4 - 1.) 657 zdnumer_so4 = aks3(ji,jj,jk) 658 zdalk_so4 = -zst * (zdnumer_so4/zdenom_so4**2) 659 660 ! HF - F : n=1, m=1 661 znumer_flu = akf3(ji,jj,jk) 662 zdenom_flu = akf3(ji,jj,jk) + zh 663 zalk_flu = zft * (znumer_flu/zdenom_flu - 1.) 664 zdnumer_flu = akf3(ji,jj,jk) 665 zdalk_flu = -zft * (zdnumer_flu/zdenom_flu**2) 666 667 ! H2O - OH 668 aphscale = 1.0 + zst/aks3(ji,jj,jk) 669 zalk_wat = akw3(ji,jj,jk)/zh - zh/aphscale 670 zdalk_wat = -akw3(ji,jj,jk)/zh**2 - 1./aphscale 671 672 ! CALCULATE [ALK]([CO3--], [HCO3-]) 673 zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil & 674 & + zalk_so4 + zalk_flu & 675 & + zalk_wat - p_alktot 676 677 zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil & 678 & + zalk_so4 + zalk_flu + zalk_wat) 679 680 zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 681 & + zdalk_so4 + zdalk_flu + zdalk_wat 682 683 ! Adapt bracketing interval 684 IF(zeqn > 0._wp) THEN 685 zh_min(ji,jj,jk) = zh_prev 686 ELSEIF(zeqn < 0._wp) THEN 687 zh_max(ji,jj,jk) = zh_prev 688 ENDIF 689 690 IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN 691 ! if the function evaluation at the current point is 692 ! not decreasing faster than with a bisection step (at least linearly) 693 ! in absolute value take one bisection step on [ph_min, ph_max] 694 ! ph_new = (ph_min + ph_max)/2d0 695 ! 696 ! In terms of [H]_new: 697 ! [H]_new = 10**(-ph_new) 698 ! = 10**(-(ph_min + ph_max)/2d0) 699 ! = SQRT(10**(-(ph_min + phmax))) 700 ! = SQRT(zh_max * zh_min) 701 zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk)) 702 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 703 ELSE 704 ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH 705 ! = -zdeqndh * LOG(10) * [H] 706 ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) 707 ! 708 ! pH_new = pH_old + \deltapH 709 ! 710 ! [H]_new = 10**(-pH_new) 711 ! = 10**(-pH_old - \Delta pH) 712 ! = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) 713 ! = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) 714 ! = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) 715 716 zh_lnfactor = -zeqn/(zdeqndh*zh_prev) 717 718 IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN 719 zh = zh_prev*EXP(zh_lnfactor) 720 ELSE 721 zh_delta = zh_lnfactor*zh_prev 722 zh = zh_prev + zh_delta 598 723 ENDIF 599 END DO 600 END DO 601 END DO 602 603 zeqn_absmin(:,:,:) = HUGE(1._wp) 604 605 DO jn = 1, jp_maxniter_atgen 606 DO jk = 1, jpk 607 DO jj = 1, jpj 608 DO ji = 1, jpi 609 IF (rmask(ji,jj,jk) == 1.) THEN 610 zfact = rhop(ji,jj,jk) / 1000. + rtrn 611 p_alktot = trb(ji,jj,jk,jptal) / zfact 612 zdic = trb(ji,jj,jk,jpdic) / zfact 613 zbot = borat(ji,jj,jk) 614 zpt = trb(ji,jj,jk,jppo4) / zfact * po4r 615 zsit = trb(ji,jj,jk,jpsil) / zfact 616 zst = sulfat (ji,jj,jk) 617 zft = fluorid(ji,jj,jk) 618 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 619 zh = zhi(ji,jj,jk) 620 zh_prev = zh 621 622 ! H2CO3 - HCO3 - CO3 : n=2, m=0 623 znumer_dic = 2.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk) 624 zdenom_dic = ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*(ak13(ji,jj,jk) + zh) 625 zalk_dic = zdic * (znumer_dic/zdenom_dic) 626 zdnumer_dic = ak13(ji,jj,jk)*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh & 627 *(4.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)) 628 zdalk_dic = -zdic*(zdnumer_dic/zdenom_dic**2) 629 630 631 ! B(OH)3 - B(OH)4 : n=1, m=0 632 znumer_bor = akb3(ji,jj,jk) 633 zdenom_bor = akb3(ji,jj,jk) + zh 634 zalk_bor = zbot * (znumer_bor/zdenom_bor) 635 zdnumer_bor = akb3(ji,jj,jk) 636 zdalk_bor = -zbot*(zdnumer_bor/zdenom_bor**2) 637 638 639 ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 640 znumer_po4 = 3.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 641 & + zh*(2.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh* ak1p3(ji,jj,jk)) 642 zdenom_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 643 & + zh*( ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh)) 644 zalk_po4 = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1 645 zdnumer_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 646 & + zh*(4.*ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 647 & + zh*(9.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 648 & + ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) & 649 & + zh*(4.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh * ak1p3(ji,jj,jk) ) ) ) 650 zdalk_po4 = -zpt * (zdnumer_po4/zdenom_po4**2) 651 652 ! H4SiO4 - H3SiO4 : n=1, m=0 653 znumer_sil = aksi3(ji,jj,jk) 654 zdenom_sil = aksi3(ji,jj,jk) + zh 655 zalk_sil = zsit * (znumer_sil/zdenom_sil) 656 zdnumer_sil = aksi3(ji,jj,jk) 657 zdalk_sil = -zsit * (zdnumer_sil/zdenom_sil**2) 658 659 ! HSO4 - SO4 : n=1, m=1 660 aphscale = 1.0 + zst/aks3(ji,jj,jk) 661 znumer_so4 = aks3(ji,jj,jk) * aphscale 662 zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh 663 zalk_so4 = zst * (znumer_so4/zdenom_so4 - 1.) 664 zdnumer_so4 = aks3(ji,jj,jk) 665 zdalk_so4 = -zst * (zdnumer_so4/zdenom_so4**2) 666 667 ! HF - F : n=1, m=1 668 znumer_flu = akf3(ji,jj,jk) 669 zdenom_flu = akf3(ji,jj,jk) + zh 670 zalk_flu = zft * (znumer_flu/zdenom_flu - 1.) 671 zdnumer_flu = akf3(ji,jj,jk) 672 zdalk_flu = -zft * (zdnumer_flu/zdenom_flu**2) 673 674 ! H2O - OH 675 aphscale = 1.0 + zst/aks3(ji,jj,jk) 676 zalk_wat = akw3(ji,jj,jk)/zh - zh/aphscale 677 zdalk_wat = -akw3(ji,jj,jk)/zh**2 - 1./aphscale 678 679 ! CALCULATE [ALK]([CO3--], [HCO3-]) 680 zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil & 681 & + zalk_so4 + zalk_flu & 682 & + zalk_wat - p_alktot 683 684 zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil & 685 & + zalk_so4 + zalk_flu + zalk_wat) 686 687 zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 688 & + zdalk_so4 + zdalk_flu + zdalk_wat 689 690 ! Adapt bracketing interval 691 IF(zeqn > 0._wp) THEN 692 zh_min(ji,jj,jk) = zh_prev 693 ELSEIF(zeqn < 0._wp) THEN 694 zh_max(ji,jj,jk) = zh_prev 695 ENDIF 696 697 IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN 698 ! if the function evaluation at the current point is 699 ! not decreasing faster than with a bisection step (at least linearly) 700 ! in absolute value take one bisection step on [ph_min, ph_max] 701 ! ph_new = (ph_min + ph_max)/2d0 702 ! 724 725 IF( zh < zh_min(ji,jj,jk) ) THEN 726 ! if [H]_new < [H]_min 727 ! i.e., if ph_new > ph_max then 728 ! take one bisection step on [ph_prev, ph_max] 729 ! ph_new = (ph_prev + ph_max)/2d0 703 730 ! In terms of [H]_new: 704 731 ! [H]_new = 10**(-ph_new) 705 ! = 10**(-(ph_min + ph_max)/2d0) 706 ! = SQRT(10**(-(ph_min + phmax))) 707 ! = SQRT(zh_max * zh_min) 708 zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk)) 709 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 710 ELSE 711 ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH 712 ! = -zdeqndh * LOG(10) * [H] 713 ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) 714 ! 715 ! pH_new = pH_old + \deltapH 716 ! 717 ! [H]_new = 10**(-pH_new) 718 ! = 10**(-pH_old - \Delta pH) 719 ! = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) 720 ! = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) 721 ! = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) 722 723 zh_lnfactor = -zeqn/(zdeqndh*zh_prev) 724 725 IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN 726 zh = zh_prev*EXP(zh_lnfactor) 727 ELSE 728 zh_delta = zh_lnfactor*zh_prev 729 zh = zh_prev + zh_delta 730 ENDIF 731 732 IF( zh < zh_min(ji,jj,jk) ) THEN 733 ! if [H]_new < [H]_min 734 ! i.e., if ph_new > ph_max then 735 ! take one bisection step on [ph_prev, ph_max] 736 ! ph_new = (ph_prev + ph_max)/2d0 737 ! In terms of [H]_new: 738 ! [H]_new = 10**(-ph_new) 739 ! = 10**(-(ph_prev + ph_max)/2d0) 740 ! = SQRT(10**(-(ph_prev + phmax))) 741 ! = SQRT([H]_old*10**(-ph_max)) 742 ! = SQRT([H]_old * zh_min) 743 zh = SQRT(zh_prev * zh_min(ji,jj,jk)) 744 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 745 ENDIF 746 747 IF( zh > zh_max(ji,jj,jk) ) THEN 748 ! if [H]_new > [H]_max 749 ! i.e., if ph_new < ph_min, then 750 ! take one bisection step on [ph_min, ph_prev] 751 ! ph_new = (ph_prev + ph_min)/2d0 752 ! In terms of [H]_new: 753 ! [H]_new = 10**(-ph_new) 754 ! = 10**(-(ph_prev + ph_min)/2d0) 755 ! = SQRT(10**(-(ph_prev + ph_min))) 756 ! = SQRT([H]_old*10**(-ph_min)) 757 ! = SQRT([H]_old * zhmax) 758 zh = SQRT(zh_prev * zh_max(ji,jj,jk)) 759 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 760 ENDIF 761 ENDIF 762 763 zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk)) 764 765 ! Stop iterations once |\delta{[H]}/[H]| < rdel 766 ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel 767 ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| 768 769 ! Alternatively: 770 ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| 771 ! ~ 1/LOG(10) * |\Delta [H]|/[H] 772 ! < 1/LOG(10) * rdel 773 774 ! Hence |zeqn/(zdeqndh*zh)| < rdel 775 776 ! rdel <-- pp_rdel_ah_target 777 l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) 778 779 IF(l_exitnow) THEN 780 rmask(ji,jj,jk) = 0. 781 ENDIF 782 783 zhi(ji,jj,jk) = zh 784 785 IF(jn >= jp_maxniter_atgen) THEN 786 zhi(ji,jj,jk) = -1._wp 787 ENDIF 788 732 ! = 10**(-(ph_prev + ph_max)/2d0) 733 ! = SQRT(10**(-(ph_prev + phmax))) 734 ! = SQRT([H]_old*10**(-ph_max)) 735 ! = SQRT([H]_old * zh_min) 736 zh = SQRT(zh_prev * zh_min(ji,jj,jk)) 737 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 789 738 ENDIF 790 END DO 791 END DO 792 END DO 739 740 IF( zh > zh_max(ji,jj,jk) ) THEN 741 ! if [H]_new > [H]_max 742 ! i.e., if ph_new < ph_min, then 743 ! take one bisection step on [ph_min, ph_prev] 744 ! ph_new = (ph_prev + ph_min)/2d0 745 ! In terms of [H]_new: 746 ! [H]_new = 10**(-ph_new) 747 ! = 10**(-(ph_prev + ph_min)/2d0) 748 ! = SQRT(10**(-(ph_prev + ph_min))) 749 ! = SQRT([H]_old*10**(-ph_min)) 750 ! = SQRT([H]_old * zhmax) 751 zh = SQRT(zh_prev * zh_max(ji,jj,jk)) 752 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 753 ENDIF 754 ENDIF 755 756 zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk)) 757 758 ! Stop iterations once |\delta{[H]}/[H]| < rdel 759 ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel 760 ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| 761 762 ! Alternatively: 763 ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| 764 ! ~ 1/LOG(10) * |\Delta [H]|/[H] 765 ! < 1/LOG(10) * rdel 766 767 ! Hence |zeqn/(zdeqndh*zh)| < rdel 768 769 ! rdel <-- pp_rdel_ah_target 770 l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) 771 772 IF(l_exitnow) THEN 773 rmask(ji,jj,jk) = 0. 774 ENDIF 775 776 zhi(ji,jj,jk) = zh 777 778 IF(jn >= jp_maxniter_atgen) THEN 779 zhi(ji,jj,jk) = -1._wp 780 ENDIF 781 782 ENDIF 783 END_3D 793 784 END DO 794 785 !
Note: See TracChangeset
for help on using the changeset viewer.
