Changeset 7162 for branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

Timestamp:

2016-11-01T14:23:51+01:00 (7 years ago)

Author:

cetlod

Message:

new top interface : Add PISCES quota model

File:

• branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90 (modified) (8 diffs)

branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

@@ -11 +11 @@
    !!                  !  2011-02  (J. Simeon, J. Orr)  Calcon salinity dependence
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Improvment of calcite dissolution
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
    !!   p4z_lys        :   Compute the CaCO3 dissolution 
@@ -18 +19 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
+   USE p4zche          !  Chemical model
    USE prtctl_trc      !  print control for debugging
    USE iom             !  I/O manager
@@ -57 +59 @@
       INTEGER, INTENT(in) ::   kt, knt ! ocean time step
       INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zalk, zdic, zph, zah2
-      REAL(wp) ::   zdispot, zfact, zcalcon, zalka, zaldi
+      REAL(wp) ::   zdispot, zfact, zcalcon
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss   
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss, zhinit, zhi, zco3sat
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_lys')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )
+      CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi, zco3sat )
       !
       zco3    (:,:,:) = 0.
       zcaldiss(:,:,:) = 0.
+      zhinit(:,:,:)   = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
       !     -------------------------------------------
       !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
       !     -------------------------------------------
-      
-      DO jn = 1, 5                               !  BEGIN OF ITERATION
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zfact = rhop(ji,jj,jk) / 1000. + rtrn
-                  zph  = hi(ji,jj,jk) * tmask(ji,jj,jk) / zfact + ( 1.-tmask(ji,jj,jk) ) * 1.e-9 ! [H+]
-                  zdic  = trb(ji,jj,jk,jpdic) / zfact
-                  zalka = trb(ji,jj,jk,jptal) / zfact
-                  ! CALCULATE [ALK]([CO3--], [HCO3-])
-                  zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph / ( aphscale(ji,jj,jk) + rtrn )  &
-                  &       + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) )
-                  ! CALCULATE [H+] and [CO3--]
-                  zaldi = zdic - zalk
-                  zah2  = SQRT( zaldi * zaldi + 4.* ( zalk * ak23(ji,jj,jk) / ak13(ji,jj,jk) ) * ( zdic + zaldi ) )
-                  zah2  = 0.5 * ak13(ji,jj,jk) / zalk * ( zaldi + zah2 )
-                  !
-                  zco3(ji,jj,jk) = zalk / ( 2. + zah2 / ak23(ji,jj,jk) ) * zfact
-                  hi(ji,jj,jk)   = zah2 * zfact
-               END DO
+
+      CALL solve_at_general(zhinit, zhi)
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
+               &                + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
+               hi(ji,jj,jk)   = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000.
             END DO
          END DO
-         !
-      END DO 
+      END DO
 
       !     ---------------------------------------------------------
@@ -111 +100 @@
                ! DEVIATION OF [CO3--] FROM SATURATION VALUE
                ! Salinity dependance in zomegaca and divide by rhop/1000 to have good units
-               zcalcon  = calcon * ( tsn(ji,jj,jk,jp_sal) / 35._wp )
+               zcalcon  = calcon * ( salinprac(ji,jj,jk) / 35._wp )
                zfact    = rhop(ji,jj,jk) / 1000._wp
                zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn )
@@ -128 +117 @@
               !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
               zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
-              zco3(ji,jj,jk)      = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk)
               !
               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk)
@@ -151 +139 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )
+      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi, zco3sat )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_lys')
@@ -170 +158 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER  ::  ji, jj, jk
       INTEGER  ::  ios                 ! Local integer output status for namelist read
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
 
       NAMELIST/nampiscal/ kdca, nca
@@ -199 +184 @@
       !
    END SUBROUTINE p4z_lys_init
-
    !!======================================================================
 END MODULE p4zlys