- Timestamp:
- 2018-12-03T12:45:01+01:00 (5 years ago)
- Location:
- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/oce_sed.F90
r10345 r10368 11 11 USE par_trc 12 12 13 USE dom_oce , ONLY : nidom => nidom !:14 13 USE dom_oce , ONLY : glamt => glamt !: longitude of t-point (degre) 15 14 USE dom_oce , ONLY : gphit => gphit !: latitude of t-point (degre) … … 21 20 USE dom_oce , ONLY : rdt => rdt !: time step for the dynamics 22 21 USE dom_oce , ONLY : nyear => nyear !: Current year 23 USE dom_oce , ONLY : nmonth => nmonth !: Current month24 USE dom_oce , ONLY : nday => nday !: Current day25 22 USE dom_oce , ONLY : ndastp => ndastp !: time step date in year/month/day aammjj 26 USE dom_oce , ONLY : nday_year => nday_year !: curent day counted from jan 1st of the current year27 23 USE dom_oce , ONLY : adatrj => adatrj !: number of elapsed days since the begining of the run 28 24 USE trc , ONLY : nittrc000 => nittrc000 … … 35 31 USE sms_pisces, ONLY : wsbio4 => wsbio4 !: sinking flux for POC 36 32 USE sms_pisces, ONLY : wsbio3 => wsbio3 !: sinking flux for GOC 37 USE sms_pisces, ONLY : wscal => wscal !: sinking flux for calcite38 33 USE sms_pisces, ONLY : wsbio2 => wsbio2 !: sinking flux for calcite 39 34 USE sms_pisces, ONLY : wsbio => wsbio !: sinking flux for calcite 35 USE sms_pisces, ONLY : ln_p5z => ln_p5z !: PISCES-QUOTA flag 40 36 USE p4zche, ONLY : akb3 => akb3 !: Chemical constants 41 37 USE sms_pisces, ONLY : ak13 => ak13 !: Chemical constants -
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedchem.F90
r10345 r10368 678 678 & - 0.07711*zsal + 0.0041249*zsal15 679 679 680 ! CONVERT FROM DIFFERENT PH SCALES 681 total2free = 1.0/(1.0 + zst/zcks) 682 free2SWS = 1. + zst/zcks + zft/(zckf*total2free) 683 total2SWS = total2free * free2SWS 684 SWS2total = 1.0 / total2SWS 685 686 680 687 ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 681 688 zak1 = 10**(zck1) * total2SWS … … 743 750 aksis(ji) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 744 751 745 ! CONVERT FROM DIFFERENT PH SCALES746 total2free = 1.0/(1.0 + zst/aks3s(ji))747 free2SWS = 1. + zst/aks3s(ji) + zft/akf3s(ji)748 total2SWS = total2free * free2SWS749 SWS2total = 1.0 / total2SWS750 751 752 ! Convert to total scale 752 753 ak1s(ji) = ak1s(ji) * SWS2total -
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/seddsr.F90
r10345 r10368 591 591 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 592 592 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 593 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 594 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn ) 593 IF ( zalpha == 0. ) THEN 594 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 595 ELSE 596 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 597 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 598 ENDIF 595 599 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 596 600 zsedtra(5) = zbeta - zsedtra(4) … … 601 605 zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) 602 606 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 603 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & 604 & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn ) 607 IF ( zalpha == 0. ) THEN 608 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 609 ELSE 610 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & 611 & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 612 ENDIF 605 613 zsedtra(1) = zalpha + 2.0 * zsedtra(2) 606 614 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) … … 609 617 zalpha = zsedtra(1) - 2.0 * zsedtra(3) 610 618 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 611 zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 612 & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn ) 619 IF ( zalpha == 0. ) THEN 620 zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0 ) 621 ELSE 622 zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 623 & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 624 ENDIF 613 625 zsedtra(1) = zalpha + 2.0 * zsedtra(3) 614 626 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) … … 617 629 zbeta = zsedtra(4) + zsedtra(6) 618 630 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 619 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & 620 & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn ) 631 IF ( zalpha == 0. ) THEN 632 zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 633 ELSE 634 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & 635 & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 636 ENDIF 621 637 zsedtra(1) = zalpha + 2.0 * zsedtra(6) 622 638 zsedtra(4) = zbeta - zsedtra(6) … … 626 642 zbeta = zsedtra(4) + zsedtra(6) 627 643 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 628 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & 629 & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn ) 644 IF ( zalpha == 0. ) THEN 645 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 646 ELSE 647 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & 648 & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 649 ENDIF 630 650 zsedtra(2) = zalpha + zsedtra(4) 631 651 zsedtra(6) = zbeta - zsedtra(4) … … 637 657 zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) 638 658 zepsi = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 639 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & 640 & + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) + rtrn ) 659 IF ( zalpha == 0. ) THEN 660 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 ) 661 ELSE 662 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & 663 & + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) ) 664 ENDIF 641 665 zsedtra(5) = zalpha + 2.0 * zsedtra(2) 642 666 zsedtra(4) = zbeta - zsedtra(5) … … 648 672 zbeta = zsedtra(4) + zsedtra(6) 649 673 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 650 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & 651 & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn ) 674 IF ( zalpha == 0. ) THEN 675 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 676 ELSE 677 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & 678 & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 679 ENDIF 652 680 zsedtra(2) = zalpha + zsedtra(4) 653 681 zsedtra(6) = zbeta - zsedtra(4) … … 657 685 zbeta = zsedtra(4) + zsedtra(6) 658 686 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 659 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & 660 & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn ) 687 IF (zalpha == 0.) THEN 688 zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. ) 689 ELSE 690 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & 691 & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 692 ENDIF 661 693 zsedtra(1) = zalpha + 2.0 * zsedtra(6) 662 694 zsedtra(4) = zbeta - zsedtra(6) … … 665 697 zalpha = zsedtra(1) - 2.0 * zsedtra(3) 666 698 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 667 zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 668 & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn ) 699 IF (zalpha == 0.) THEN 700 zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0) 701 ELSE 702 zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 703 & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 704 ENDIF 669 705 zsedtra(1) = zalpha + 2.0 * zsedtra(3) 670 706 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) … … 673 709 zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) 674 710 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 675 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & 676 & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn ) 711 IF ( zalpha == 0. ) THEN 712 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 713 ELSE 714 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & 715 & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 716 ENDIF 677 717 zsedtra(1) = zalpha + 2.0 * zsedtra(2) 678 718 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) … … 683 723 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 684 724 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 685 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 686 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn ) 725 IF ( zalpha == 0. ) THEN 726 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 727 ELSE 728 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 729 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 730 ENDIF 687 731 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 688 732 zsedtra(5) = zbeta - zsedtra(4) -
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/seddta.F90
r10345 r10368 49 49 50 50 REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag 51 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3 , zwscal51 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3 52 52 REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep 53 53 … … 97 97 zwsbio4(ji,jj) = wsbio2 / rday 98 98 zwsbio3(ji,jj) = wsbio / rday 99 zwscal (ji,jj) = wsbio2 / rday100 99 END DO 101 100 END DO … … 106 105 zdep = e3t_n(ji,jj,ikt) / r2dttrc 107 106 zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) / rday ) 108 zwscal (ji,jj) = MIN( 0.99 * zdep, wscal (ji,jj,ikt) / rday )109 107 zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) / rday ) 110 108 END DO … … 130 128 trc_data(ji,jj,12 ) = MIN(trb(ji,jj,ikt,jppoc), 1E-4) * zwsbio3(ji,jj) * 1E3 131 129 trc_data(ji,jj,13 ) = MIN(trb(ji,jj,ikt,jpgoc), 1E-4) * zwsbio4(ji,jj) * 1E3 132 trc_data(ji,jj,14) = MIN(trb(ji,jj,ikt,jpcal), 1E-4) * zws cal(ji,jj) * 1E3130 trc_data(ji,jj,14) = MIN(trb(ji,jj,ikt,jpcal), 1E-4) * zwsbio4(ji,jj) * 1E3 133 131 trc_data(ji,jj,15) = tsn(ji,jj,ikt,jp_tem) 134 132 trc_data(ji,jj,16) = tsn(ji,jj,ikt,jp_sal) -
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedini.F90
r10345 r10368 472 472 ENDIF 473 473 474 IF ( ln_p5z .AND. ln_sed_2way ) CALL ctl_stop( '2 ways coupling with sediment cannot be activated with PISCES-QUOTA' ) 475 474 476 REWIND( numnamsed_ref ) ! Namelist nam_geom in reference namelist : Pisces variables 475 477 READ ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903)
Note: See TracChangeset
for help on using the changeset viewer.