New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6971 – NEMO

Changeset 6971


Ignore:
Timestamp:
2016-10-03T09:52:43+02:00 (8 years ago)
Author:
clem
Message:

update from 3.6 stable

Location:
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

    r6861 r6971  
    26822682      !!---------------------------------------------------------------------- 
    26832683      ! 
    2684       ALLOCATE( ztab(jpiglo,4,num_fields), znorthloc(jpi,4,num_fields), zfoldwk(jpi,4,num_fields), znorthgloio(jpi,4,num_fields,jpni) )   ! expanded to 3 dimensions 
     2684      ALLOCATE( ztab(jpiglo,4,num_fields), znorthloc(jpi,4,num_fields), zfoldwk(jpi,4,num_fields),   &  
     2685            &   znorthgloio(jpi,4,num_fields,jpni) )   ! expanded to 3 dimensions 
    26852686      ALLOCATE( ztabl(jpi,4,num_fields), ztabr(jpi*jpmaxngh, 4,num_fields) ) 
    26862687      ! 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/step.F90

    r6970 r6971  
    353353      ENDIF 
    354354#endif 
    355       IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
     355      IF( lk_diaobs        )   CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
    356356 
    357357      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90

    r6287 r6971  
    3131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sio3eq   ! chemistry of Si 
    3232   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fekeq    ! chemistry of Fe 
    33    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ::   chemc    ! Solubilities of O2 and CO2 
     33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemc    ! Solubilities of O2 and CO2 
    3434   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2   ! Solubilities of O2 and CO2 
     35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature 
    3536 
    3637   REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm 
     
    3940   REAL(wp) ::   o2atm  = 1. / ( 1000. * 0.20946 )   
    4041 
    41    REAL(wp) ::   akcc1  = -171.9065       ! coeff. for apparent solubility equilibrium 
    42    REAL(wp) ::   akcc2  =   -0.077993     ! Millero et al. 1995 from Mucci 1983 
    43    REAL(wp) ::   akcc3  = 2839.319         
    44    REAL(wp) ::   akcc4  =   71.595         
    45    REAL(wp) ::   akcc5  =   -0.77712       
    46    REAL(wp) ::   akcc6  =    0.00284263    
    47    REAL(wp) ::   akcc7  =  178.34         
    48    REAL(wp) ::   akcc8  =   -0.07711      
    49    REAL(wp) ::   akcc9  =    0.0041249    
    50  
    51    REAL(wp) ::   rgas   = 83.143         ! universal gas constants 
     42   REAL(wp) ::   rgas   = 83.14472       ! universal gas constants 
    5243   REAL(wp) ::   oxyco  = 1. / 22.4144   ! converts from liters of an ideal gas to moles 
    5344 
    5445   REAL(wp) ::   bor1   = 0.00023        ! borat constants 
    5546   REAL(wp) ::   bor2   = 1. / 10.82 
    56  
    57    REAL(wp) ::   ca0    = -162.8301      ! WEISS & PRICE 1980, units mol/(kg atm) 
    58    REAL(wp) ::   ca1    =  218.2968 
    59    REAL(wp) ::   ca2    =   90.9241 
    60    REAL(wp) ::   ca3    =   -1.47696 
    61    REAL(wp) ::   ca4    =    0.025695 
    62    REAL(wp) ::   ca5    =   -0.025225 
    63    REAL(wp) ::   ca6    =    0.0049867 
    64  
    65    REAL(wp) ::   c10    = -3670.7        ! Coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)    
    66    REAL(wp) ::   c11    =    62.008      
    67    REAL(wp) ::   c12    =    -9.7944     
    68    REAL(wp) ::   c13    =     0.0118      
    69    REAL(wp) ::   c14    =    -0.000116 
    70  
    71    REAL(wp) ::   c20    = -1394.7       ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)    
    72    REAL(wp) ::   c21    =    -4.777    
    73    REAL(wp) ::   c22    =     0.0184    
    74    REAL(wp) ::   c23    =    -0.000118 
    7547 
    7648   REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate 
     
    146118      REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 
    147119      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2 
    148       REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1 
    149       REAL(wp) ::   zis  , zis2 , zsal15, zisqrt 
     120      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat 
     121      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1  , za2 
    150122      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw 
    151123      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1 
     
    154126      IF( nn_timing == 1 )  CALL timing_start('p4z_che') 
    155127      ! 
     128      ! Computations of chemical constants require in situ temperature 
     129      ! Here a quite simple formulation is used to convert  
     130      ! potential temperature to in situ temperature. The errors is less than  
     131      ! 0.04°C relative to an exact computation 
     132      ! --------------------------------------------------------------------- 
     133      DO jk = 1, jpk 
     134         DO jj = 1, jpj 
     135            DO ji = 1, jpi 
     136               zpres = fsdept(ji,jj,jk) / 1000. 
     137               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (tsn(ji,jj,jk,jp_sal) - 35.0) ) 
     138               za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 ) 
     139               tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2 
     140            END DO 
     141         END DO 
     142      END DO 
     143      ! 
    156144      ! CHEMICAL CONSTANTS - SURFACE LAYER 
    157145      ! ---------------------------------- 
     
    161149         DO ji = 1, jpi 
    162150            !                             ! SET ABSOLUTE TEMPERATURE 
    163             ztkel = tsn(ji,jj,1,jp_tem) + 273.15 
     151            ztkel = tempis(ji,jj,1) + 273.15 
    164152            zt    = ztkel * 0.01 
    165153            zt2   = zt * zt 
     
    169157            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 
    170158            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
    171             zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 ) 
     159            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    & 
     160            &       + 0.0047036e-4*ztkel**2) 
    172161            !                             ! SET SOLUBILITIES OF O2 AND CO2  
    173             chemc(ji,jj) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm) 
     162            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm) 
     163            chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 
     164            chemc(ji,jj,3) = 57.7 - 0.118*ztkel 
    174165            ! 
    175166         END DO 
     
    184175!CDIR NOVERRCHK 
    185176            DO ji = 1, jpi 
    186               ztkel = tsn(ji,jj,jk,jp_tem) + 273.15 
     177              ztkel = tempis(ji,jj,jk) + 273.15 
    187178              zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35. 
    188179              zsal2 = zsal * zsal 
    189               ztgg  = LOG( ( 298.15 - tsn(ji,jj,jk,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
     180              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
    190181              ztgg2 = ztgg  * ztgg 
    191182              ztgg3 = ztgg2 * ztgg 
     
    210201            DO ji = 1, jpi 
    211202 
    212                ! SET PRESSION 
    213                zpres   = 1.025e-1 * fsdept(ji,jj,jk) 
     203               ! SET PRESSION ACCORDING TO SAUNDER (1980) 
     204               zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) ) 
     205               zc1 = 5.92E-3 + zplat**2 * 5.25E-3 
     206               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*fsdept(ji,jj,jk)))) / 4.42E-6 
     207               zpres = zpres / 10.0 
    214208 
    215209               ! SET ABSOLUTE TEMPERATURE 
    216                ztkel   = tsn(ji,jj,jk,jp_tem) + 273.15 
     210               ztkel   = tempis(ji,jj,jk) + 273.15 
    217211               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35. 
    218212               zsqrt  = SQRT( zsal ) 
     
    223217               zis2   = zis * zis 
    224218               zisqrt = SQRT( zis ) 
    225                ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20. 
     219               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20. 
    226220 
    227221               ! CHLORINITY (WOOSTER ET AL., 1969) 
     
    256250 
    257251 
    258                zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal 
    259                zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2 
     252               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO  
     253               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 
     254               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  & 
     255                  - 0.011555*zsal + 0.0001152*zsal*zsal) 
     256               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      & 
     257                  - 0.01781*zsal + 0.0001122*zsal*zsal) 
    260258 
    261259               ! PKW (H2O) (DICKSON AND RILEY, 1979) 
     
    266264               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
    267265               !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 
    268                zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   & 
    269                   &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15 
     266               zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   & 
     267                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  & 
     268                  &      - 0.07711*zsal + 0.0041249*zsal15 
    270269 
    271270               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 
     
    337336      !!                     ***  ROUTINE p4z_che_alloc  *** 
    338337      !!---------------------------------------------------------------------- 
    339       ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj), chemo2(jpi,jpj,jpk),   & 
    340       &         STAT=p4z_che_alloc ) 
     338      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk),   & 
     339      &         tempis(jpi,jpj,jpk), STAT=p4z_che_alloc ) 
    341340      ! 
    342341      IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90

    r6287 r6971  
    8686      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan 
    8787      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact 
     88      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff 
    8889      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 
    8990      REAL(wp) ::   zyr_dec, zdco2dt 
    9091      CHARACTER (len=25) :: charout 
    91       REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d  
     92      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d, zpco2atm  
    9293      !!--------------------------------------------------------------------- 
    9394      ! 
    9495      IF( nn_timing == 1 )  CALL timing_start('p4z_flx') 
    9596      ! 
    96       CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx ) 
     97      CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 
    9798      ! 
    9899 
     
    183184      DO jj = 1, jpj 
    184185         DO ji = 1, jpi 
     186            ztkel  = tsn(ji,jj,1,jp_tem) + 273.15 
     187            zsal   = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35. 
     188            zvapsw = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal) 
     189            zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) 
     190            zxc2 = (1.0 - zpco2atm(ji,jj) * 1E-6 )**2 
     191            zfugcoeff = EXP(patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   & 
     192            &           / (82.05736 * ztkel)) 
     193            zfco2 = zpco2atm(ji,jj) * zfugcoeff 
     194 
    185195            ! Compute CO2 flux for the sea and air 
    186             zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj) * zkgco2(ji,jj)   ! (mol/L) * (m/s) 
    187             zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ? 
     196            zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s) 
     197            zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ? 
    188198            oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 
    189199            ! compute the trend 
    190             tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1) 
     200            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1) * tmask(ji,jj,1) 
    191201 
    192202            ! Compute O2 flux  
    193             zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
    194             zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj) 
    195             zoflx(ji,jj) = zfld16 - zflu16 
     203            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
     204            zflu16 = trb(ji,jj,1,jpoxy) * zkgo2(ji,jj) 
     205            zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1) 
    196206            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / fse3t(ji,jj,1) 
    197207         END DO 
     
    224234         ENDIF 
    225235         IF( iom_use( "Dpco2" ) ) THEN 
    226            zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1) 
     236           zw2d(:,:) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
    227237           CALL iom_put( "Dpco2" ,  zw2d ) 
    228238         ENDIF 
     
    240250            trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1)  
    241251            trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1)  
    242             trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1)  
    243          ENDIF 
    244       ENDIF 
    245       ! 
    246       CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx ) 
     252            trc2d(:,:,jp_pcs0_2d + 3) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
     253         ENDIF 
     254      ENDIF 
     255      ! 
     256      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 
    247257      ! 
    248258      IF( nn_timing == 1 )  CALL timing_stop('p4z_flx') 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90

    r6204 r6971  
    4444   REAL(wp), PUBLIC ::  xkdoc       !:  2nd half-sat. of DOC remineralization   
    4545   REAL(wp), PUBLIC ::  concbfe     !:  Fe half saturation for bacteria  
     46   REAL(wp), PUBLIC ::  oxymin      !:  half saturation constant for anoxia 
    4647   REAL(wp), PUBLIC ::  qnfelim     !:  optimal Fe quota for nanophyto 
    4748   REAL(wp), PUBLIC ::  qdfelim     !:  optimal Fe quota for diatoms 
     
    121122               zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) 
    122123               zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concbnh4 ) 
    123                zlim3    = trb(ji,jj,jk,jpfer) / ( concbfe + trb(ji,jj,jk,jpfer) ) 
     124               zlim3    = biron(ji,jj,jk)     / ( concbfe + biron(ji,jj,jk) ) 
    124125               zlim4    = trb(ji,jj,jk,jpdoc) / ( xkdoc   + trb(ji,jj,jk,jpdoc) ) 
    125126               xlimbacl(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 ) 
     
    187188      END DO 
    188189      ! 
     190      DO jk = 1, jpkm1 
     191         DO jj = 1, jpj 
     192            DO ji = 1, jpi 
     193               ! denitrification factor computed from O2 levels 
     194               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    & 
     195                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  ) 
     196               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) ) 
     197            END DO 
     198         END DO 
     199      END DO 
    189200      ! 
    190201      IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics 
     
    216227      NAMELIST/nampislim/ concnno3, concdno3, concnnh4, concdnh4, concnfer, concdfer, concbfe,   & 
    217228         &                concbno3, concbnh4, xsizedia, xsizephy, xsizern, xsizerd,          &  
    218          &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r 
     229         &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r, oxymin 
    219230      INTEGER :: ios                 ! Local integer output status for namelist read 
    220231 
     
    249260         WRITE(numout,*) '    Minimum size criteria for nanophyto      xsizephy  = ', xsizephy 
    250261         WRITE(numout,*) '    Fe half saturation for bacteria          concbfe   = ', concbfe 
     262         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin   =' , oxymin 
    251263         WRITE(numout,*) '    optimal Fe quota for nano.               qnfelim   = ', qnfelim 
    252264         WRITE(numout,*) '    Optimal Fe quota for diatoms             qdfelim   = ', qdfelim 
    253265      ENDIF 
    254  
     266      ! 
     267      nitrfac (:,:,:) = 0._wp 
     268      ! 
    255269   END SUBROUTINE p4z_lim_init 
    256270 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

    r6287 r6971  
    6565      REAL(wp) ::   zomegaca, zexcess, zexcess0 
    6666      CHARACTER (len=25) :: charout 
    67       REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss    
     67      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss    
    6868      !!--------------------------------------------------------------------- 
    6969      ! 
    7070      IF( nn_timing == 1 )  CALL timing_start('p4z_lys') 
    7171      ! 
    72       CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss ) 
     72      CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss ) 
    7373      ! 
    7474      zco3    (:,:,:) = 0. 
     
    120120               zcalcon  = calcon * ( tsn(ji,jj,jk,jp_sal) / 35._wp ) 
    121121               zfact    = rhop(ji,jj,jk) / 1000._wp 
    122                zomegaca = ( zcalcon * zco3(ji,jj,jk) * zfact ) / aksp(ji,jj,jk)  
     122               zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn ) 
     123               zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn ) 
    123124 
    124125               ! SET DEGREE OF UNDER-/SUPERSATURATION 
     
    149150      IF( lk_iomput .AND. knt == nrdttrc ) THEN 
    150151         IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )          * tmask(:,:,:) ) 
    151          IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:) * 1.e+3               * tmask(:,:,:) ) 
    152          IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", aksp(:,:,:) * 1.e+3 / calcon      * tmask(:,:,:) ) 
    153          IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) ) 
     152         IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:)    * 1.e+3            * tmask(:,:,:) ) 
     153         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", zco3sat(:,:,:) * 1.e+3            * tmask(:,:,:) ) 
     154         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) ) 
    154155      ELSE 
    155156         IF( ln_diatrc ) THEN 
    156157            trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:) 
    157158            trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:) 
    158             trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:) 
     159            trc3d(:,:,:,jp_pcs0_3d + 2) = zco3sat(:,:,:)           * tmask(:,:,:) 
    159160         ENDIF 
    160161      ENDIF 
     
    166167      ENDIF 
    167168      ! 
    168       CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss ) 
     169      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss ) 
    169170      ! 
    170171      IF( nn_timing == 1 )  CALL timing_stop('p4z_lys') 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90

    r6746 r6971  
    114114      !                                        !  -------------------------------------- 
    115115      IF( l_trcdm2dc ) THEN                     !  diurnal cycle 
    116          !                                       ! 1% of qsr to compute euphotic layer 
    117          zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr 
    118116         ! 
    119117         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
    120118         ! 
    121          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )  
     119         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )  
    122120         ! 
    123121         DO jk = 1, nksrp       
     
    136134         ! 
    137135      ELSE 
    138          ! 1% of qsr to compute euphotic layer 
    139          zqsr100(:,:) = 0.01 * qsr(:,:) 
    140136         ! 
    141137         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
    142138         ! 
    143          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )  
     139         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )  
    144140         ! 
    145141         DO jk = 1, nksrp       
     
    169165         DO jj = 1, jpj 
    170166           DO ji = 1, jpi 
    171               IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN 
     167              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= zqsr100(ji,jj) )  THEN 
    172168                 neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer 
    173169                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint 
     
    242238   END SUBROUTINE p4z_opt 
    243239 
    244    SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 )  
     240   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0, pqsr100 )  
    245241      !!---------------------------------------------------------------------- 
    246242      !!                  ***  routine p4z_opt_par  *** 
     
    251247      !!---------------------------------------------------------------------- 
    252248      !! * arguments 
    253       INTEGER, INTENT(in)                                       ::  kt            !   ocean time-step 
    254       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)              ::  pqsr          !   shortwave 
    255       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B) 
    256       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0   
     249      INTEGER, INTENT(in)                                        ::  kt            !   ocean time-step 
     250      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)               ::  pqsr          !   shortwave 
     251      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)            ::  pe1 , pe2 , pe3   !  PAR ( R-G-B) 
     252      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL  ::  pe0  
     253      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(out)  , OPTIONAL  ::  pqsr100   
    257254      !! * local variables 
    258255      INTEGER    ::   ji, jj, jk     ! dummy loop indices 
     
    264261      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:) 
    265262      ENDIF 
     263 
     264      !  Light at the euphotic depth  
     265      IF( PRESENT( pqsr100 ) )  pqsr100(:,:) = 0.01 * 3. * zqsr(:,:) 
    266266      ! 
    267267      IF( PRESENT( pe0 ) ) THEN     !  W-level 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90

    r6746 r6971  
    136136                  zval = MAX( 1., zstrn(ji,jj) ) 
    137137                  zval = 1.5 * zval / ( 12. + zval ) 
    138                   zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval * ( 1. - fr_i(ji,jj) ) 
     138                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval 
    139139                  zprdia(ji,jj,jk) = zprbio(ji,jj,jk) 
    140140               ENDIF 
     
    202202                      zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj) 
    203203                      ! 
    204                       zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -znanotot ) ) 
    205                       zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trb(ji,jj,jk,jpdia) + rtrn ) 
    206  
    207                       zpislopen =  zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch)                & 
    208                         &          / ( trb(ji,jj,jk,jpphy) * 12.                  + rtrn )   & 
    209                         &          / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn ) 
    210  
    211                       zpislope2n = zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch)                & 
    212                         &          / ( trb(ji,jj,jk,jpdia) * 12.                  + rtrn )   & 
    213                         &          / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn ) 
     204                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -znanotot ) )           & 
     205                         &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn) 
     206                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trb(ji,jj,jk,jpdia) + rtrn )   & 
     207                         &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn) 
    214208 
    215209                      ! Computation of production function for Carbon 
    216210                      !  --------------------------------------------- 
     211                      zpislopen  =  zpislopead(ji,jj,jk)  / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn ) 
     212                      zpislope2n =  zpislopead2(ji,jj,jk) / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn ) 
    217213                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) ) 
    218214                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) ) 
     
    220216                      !  Computation of production function for Chlorophyll 
    221217                      !-------------------------------------------------- 
    222                       zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) ) 
    223                       zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) ) 
     218                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) ) 
     219                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) ) 
    224220                  ENDIF 
    225221               END DO 
     
    227223         END DO 
    228224      ENDIF 
    229  
    230  
     225       
    231226      !  Computation of a proxy of the N/C ratio 
    232227      !  --------------------------------------- 
     
    278273            zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) ) 
    279274            zmxlday = zmxltst * zmxltst * r1_rday 
    280             zmixnano(ji,jj) = 1. - zmxlday / ( 2. + zmxlday ) 
    281             zmixdiat(ji,jj) = 1. - zmxlday / ( 4. + zmxlday ) 
     275            zmixnano(ji,jj) = 1. - zmxlday / ( 1. + zmxlday ) 
     276            zmixdiat(ji,jj) = 1. - zmxlday / ( 2. + zmxlday ) 
    282277         END DO 
    283278      END DO 
    284279  
    285       !  Mixed-layer effect on production                                                                                
     280      !  Mixed-layer effect on production  
     281      !  Sea-ice effect on production 
     282 
    286283      DO jk = 1, jpkm1 
    287284         DO jj = 1, jpj 
     
    291288                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj) 
    292289               ENDIF 
     290                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
     291                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
    293292            END DO 
    294293         END DO 
     
    330329      END DO 
    331330 
    332       IF( ln_newprod ) THEN 
    333 !CDIR NOVERRCHK 
    334          DO jk = 1, jpkm1 
    335 !CDIR NOVERRCHK 
    336             DO jj = 1, jpj 
    337 !CDIR NOVERRCHK 
    338                DO ji = 1, jpi 
    339                   IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
    340                      zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj) 
    341                      zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj) 
    342                   ENDIF 
    343                   IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    344                      !  production terms for nanophyto. ( chlorophyll ) 
    345                      znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 
    346                      zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 
    347                      zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
    348                      zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / & 
    349                                         & (  zpislopead(ji,jj,jk) * znanotot +rtrn) 
    350                      !  production terms for diatomees ( chlorophyll ) 
    351                      zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 
    352                      zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 
    353                      zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
    354                      zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / & 
    355                                         & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn ) 
    356                   ENDIF 
    357                END DO 
    358             END DO 
    359          END DO 
    360       ELSE 
    361 !CDIR NOVERRCHK 
    362          DO jk = 1, jpkm1 
    363 !CDIR NOVERRCHK 
    364             DO jj = 1, jpj 
    365 !CDIR NOVERRCHK 
    366                DO ji = 1, jpi 
    367                   IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    368                      !  production terms for nanophyto. ( chlorophyll ) 
    369                      znanotot = enano(ji,jj,jk) 
    370                      zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trb(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk) 
    371                      zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
    372                      zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 144. * zprod            & 
    373                      &                    / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * znanotot +rtrn ) 
    374                      !  production terms for diatomees ( chlorophyll ) 
    375                      zdiattot = ediat(ji,jj,jk) 
    376                      zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trb(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk) 
    377                      zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
    378                      zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 144. * zprod             & 
    379                      &                    / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn ) 
    380                   ENDIF 
    381                END DO 
    382             END DO 
    383          END DO 
    384       ENDIF 
     331!CDIR NOVERRCHK 
     332      DO jk = 1, jpkm1 
     333!CDIR NOVERRCHK 
     334         DO jj = 1, jpj 
     335!CDIR NOVERRCHK 
     336            DO ji = 1, jpi 
     337               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
     338                  zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj) 
     339                  zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj) 
     340               ENDIF 
     341               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
     342                  !  production terms for nanophyto. ( chlorophyll ) 
     343                  znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 
     344                  zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 
     345                  zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
     346                  zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / & 
     347                                     & (  zpislopead(ji,jj,jk) * znanotot +rtrn) 
     348                  !  production terms for diatomees ( chlorophyll ) 
     349                  zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 
     350                  zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 
     351                  zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
     352                  zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / & 
     353                                     & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn ) 
     354               ENDIF 
     355            END DO 
     356         END DO 
     357      END DO 
    385358 
    386359      !   Update the arrays TRA which contain the biological sources and sinks 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90

    r5385 r6971  
    4444   REAL(wp), PUBLIC ::  xsiremlab  !: fast remineralisation rate of POC  
    4545   REAL(wp), PUBLIC ::  xsilab     !: fraction of labile biogenic silica  
    46    REAL(wp), PUBLIC ::  oxymin     !: halk saturation constant for anoxia  
    47  
    4846 
    4947   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitr     !: denitrification array 
     
    111109                  zdepprod(ji,jj,jk) = zdepmin**0.273 
    112110               ENDIF 
    113             END DO 
    114          END DO 
    115       END DO 
    116  
    117       DO jk = 1, jpkm1 
    118          DO jj = 1, jpj 
    119             DO ji = 1, jpi 
    120                ! denitrification factor computed from O2 levels 
    121                nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    & 
    122                   &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  ) 
    123                nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) ) 
    124111            END DO 
    125112         END DO 
     
    357344      !! 
    358345      !!---------------------------------------------------------------------- 
    359       NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab,   & 
    360       &                   oxymin 
     346      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab 
    361347      INTEGER :: ios                 ! Local integer output status for namelist read 
    362348 
     
    380366         WRITE(numout,*) '    fraction of labile biogenic silica        xsilab    =', xsilab 
    381367         WRITE(numout,*) '    NH4 nitrification rate                    nitrif    =', nitrif 
    382          WRITE(numout,*) '    halk saturation constant for anoxia       oxymin    =', oxymin 
    383368      ENDIF 
    384369      ! 
    385       nitrfac (:,:,:) = 0._wp 
    386370      denitr  (:,:,:) = 0._wp 
    387371      denitnh4(:,:,:) = 0._wp 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90

    r6204 r6971  
    159159      IF( ln_ndepo ) THEN 
    160160         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN 
    161             CALL fld_read( kt, 1, sf_ndepo ) 
    162             DO jj = 1, jpj 
    163                DO ji = 1, jpi 
    164                   nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / rno3 / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn ) 
    165                END DO 
    166             END DO 
     161             zcoef = rno3 * 14E6 * ryyss 
     162             CALL fld_read( kt, 1, sf_ndepo ) 
     163             nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / fse3t(:,:,1)  
     164         ENDIF 
     165         IF( lk_vvl ) THEN 
     166           zcoef = rno3 * 14E6 * ryyss 
     167           nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / fse3t(:,:,1)  
    167168         ENDIF 
    168169      ENDIF 
     
    455456            DO jj = 1, jpj 
    456457               DO ji = 1, jpi 
    457                   zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) ) 
     458                  zexpide   = MIN( 8.,( gdept_0(ji,jj,jk) / 500. )**(-1.5) ) 
    458459                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2 
    459460                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 ) 
     
    465466         ironsed(:,:,jpk) = 0._wp 
    466467         DO jk = 1, jpkm1 
    467             ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday ) 
     468            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday ) 
    468469         END DO 
    469470         DEALLOCATE( zcmask) 
     
    483484         CALL iom_close( numhydro ) 
    484485         ! 
    485          hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp 
     486         DO jk = 1, jpk 
     487            hydrofe(:,:,jk) = ( hydrofe(:,:,jk) * hratio ) / ( e1e2t(:,:) * e3t_0(:,:,jk) * ryyss + rtrn ) / 1000._wp 
     488         ENDDO 
    486489         ! 
    487490      ENDIF 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90

    r6308 r6971  
    102102         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    103103 
    104          IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file 
     104         IF( ln_rsttr .AND. .NOT.ln_top_euler .AND.   &                     ! Restart: read in restart  file 
    105105            iom_varid( numrtr, 'sbc_'//TRIM(ctrcnm(1))//'_b', ldstop = .FALSE. ) > 0 ) THEN 
    106106            IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc surface tracer content forcing fields red in the restart file' 
     
    190190      !                                           Write in the tracer restar  file 
    191191      !                                          ******************************* 
    192       IF( lrst_trc ) THEN 
     192      IF( lrst_trc .AND. .NOT.ln_top_euler ) THEN 
    193193         IF(lwp) WRITE(numout,*) 
    194194         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in tracer restart file ',   & 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90

    r6308 r6971  
    6868         IF( ln_trcdmp )        CALL trc_dmp( kstp )            ! internal damping trends 
    6969                                CALL trc_adv( kstp )            ! horizontal & vertical advection  
     70         IF( ln_zps ) THEN 
     71           IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kstp, jptra, trb, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom 
     72           ELSE                 ; CALL zps_hde    ( kstp, jptra, trb, gtru, gtrv )                                      !  only bottom 
     73           ENDIF 
     74         ENDIF 
    7075                                CALL trc_ldf( kstp )            ! lateral mixing 
    7176         IF( .NOT. lk_offline .AND. lk_zdfkpp )    & 
     
    7580#endif 
    7681                                CALL trc_zdf( kstp )            ! vertical mixing and after tracer fields 
     82         ! 
    7783                                CALL trc_nxt( kstp )            ! tracer fields at next time step      
    7884         IF( ln_trcrad )        CALL trc_rad( kstp )            ! Correct artificial negative concentrations 
     
    8389#endif 
    8490 
    85          IF( ln_zps  .AND. .NOT. ln_isfcav)        & 
    86             &            CALL zps_hde    ( kstp, jptra, trn, gtru, gtrv )   ! Partial steps: now horizontal gradient of passive 
    87          IF( ln_zps .AND.        ln_isfcav)        & 
    88             &            CALL zps_hde_isf( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! Partial steps: now horizontal gradient of passive 
    89                                                                 ! tracers at the bottom ocean level 
    90          ! 
    9191      ELSE                                               ! 1D vertical configuration 
    9292                                CALL trc_sbc( kstp )            ! surface boundary condition 
     
    100100      ! 
    101101      IF( nn_timing == 1 )   CALL timing_stop('trc_trp') 
     102      ! 
     1039400  FORMAT(a25,i4,D23.16) 
    102104      ! 
    103105   END SUBROUTINE trc_trp 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r6746 r6971  
    2626   USE trcdta          ! initialisation from files 
    2727   USE daymod          ! calendar manager 
    28    USE zpshde          ! partial step: hor. derivative   (zps_hde routine) 
    2928   USE prtctl_trc      ! Print control passive tracers (prt_ctl_trc_init routine) 
    3029   USE trcsub          ! variables to substep passive tracers 
     
    142141  
    143142      tra(:,:,:,:) = 0._wp 
    144       IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav )   &              ! Partial steps: before horizontal gradient of passive 
    145         &    CALL zps_hde    ( nit000, jptra, trn, gtru, gtrv  )  ! Partial steps: before horizontal gradient 
    146       IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav )   & 
    147         &    CALL zps_hde_isf( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )       ! tracers at the bottom ocean level 
    148  
    149  
    150143      ! 
    151144      IF( nn_dttrc /= 1 )        CALL trc_sub_ini      ! Initialize variables for substepping passive tracers 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcstp.F90

    r6204 r6971  
    3333   REAL(wp) :: rdt_sampl 
    3434   INTEGER  :: nb_rec_per_day 
    35    INTEGER  :: isecfst, iseclast 
     35   REAL(wp) :: rsecfst, rseclast 
    3636   LOGICAL  :: llnew 
    3737 
     
    5959      REAL(wp)              ::  ztrai 
    6060      CHARACTER (len=25)    ::  charout  
    61  
    6261      !!------------------------------------------------------------------- 
    6362      ! 
     
    9493                                   CALL trc_sms      ( kt )       ! tracers: sinks and sources 
    9594                                   CALL trc_trp      ( kt )       ! transport of passive tracers 
     95 
    9696         IF( kt == nittrc000 ) THEN 
    9797            CALL iom_close( numrtr )       ! close input tracer restart file 
     
    105105      ENDIF 
    106106      ! 
     107 
    107108      ztrai = 0._wp                                                   !  content of all tracers 
    108109      DO jn = 1, jptra 
     
    110111      END DO 
    111112      IF( lwp ) WRITE(numstr,9300) kt,  ztrai / areatot 
    112 9300  FORMAT(i10,e18.10) 
     1139300  FORMAT(i10,D23.16) 
    113114      ! 
    114115      IF( nn_timing == 1 )   CALL timing_stop('trc_stp') 
     
    130131      INTEGER, INTENT(in) ::   kt 
    131132      INTEGER  :: jn 
     133      REAL(wp) :: zkt 
     134      CHARACTER(len=1)               ::   cl1                      ! 1 character 
     135      CHARACTER(len=2)               ::   cl2                      ! 2 characters 
    132136 
    133137      IF( kt == nittrc000 ) THEN 
    134138         IF( ln_cpl )  THEN   
    135             rdt_sampl = 86400. / ncpl_qsr_freq 
     139            rdt_sampl = rday / ncpl_qsr_freq 
    136140            nb_rec_per_day = ncpl_qsr_freq 
    137141         ELSE   
    138             rdt_sampl = MAX( 3600., rdt * nn_dttrc ) 
    139             nb_rec_per_day = INT( 86400 / rdt_sampl ) 
     142            rdt_sampl = MAX( 3600., rdttrc(1) ) 
     143            nb_rec_per_day = INT( rday / rdt_sampl ) 
    140144         ENDIF 
    141145         ! 
     
    146150         ENDIF 
    147151         ! 
     152         ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_day ) ) 
     153         ! 
    148154         !                                            !* Restart: read in restart file 
    149          IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean', ldstop = .FALSE. ) > 0 ) THEN  
    150             IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file' 
     155         IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 .AND. & 
     156                            iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 .AND. & 
     157                            iom_varid( numrtr, 'ktdcy'    , ldstop = .FALSE. ) > 0 ) THEN  
     158            CALL iom_get( numrtr, 'ktdcy', zkt )   !  A mean of qsr 
     159            rsecfst = INT( zkt ) * rdttrc(1) 
     160            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s ' 
    151161            CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean )   !  A mean of qsr 
     162            DO jn = 1, nb_rec_per_day  
     163             IF( jn <= 9 )  THEN 
     164               WRITE(cl1,'(i1)') jn 
     165               CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr 
     166             ELSE 
     167               WRITE(cl2,'(i2.2)') jn 
     168               CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr 
     169             ENDIF 
     170           ENDDO 
    152171         ELSE                                         !* no restart: set from nit000 values 
    153172            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean set to nit000 values' 
     173            rsecfst  = kt * rdttrc(1) 
     174            ! 
    154175            qsr_mean(:,:) = qsr(:,:) 
    155          ENDIF 
    156          ! 
    157          ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_day ) ) 
    158          DO jn = 1, nb_rec_per_day 
    159              qsr_arr(:,:,jn) = qsr_mean(:,:) 
    160          ENDDO 
    161          ! 
    162          isecfst  = nsec_year + nsec1jan000   !   number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 
    163          iseclast = isecfst 
    164          ! 
    165       ENDIF 
    166       ! 
    167       iseclast = nsec_year + nsec1jan000 
    168       llnew   = ( iseclast - isecfst )  > INT( rdt_sampl )   !   new shortwave to store 
    169       IF( kt /= nittrc000 .AND. llnew ) THEN 
     176            DO jn = 1, nb_rec_per_day 
     177               qsr_arr(:,:,jn) = qsr_mean(:,:) 
     178            ENDDO 
     179         ENDIF 
     180         ! 
     181      ENDIF 
     182      ! 
     183      rseclast = kt * rdttrc(1) 
     184      ! 
     185      llnew   = ( rseclast - rsecfst ) .ge.  rdt_sampl    !   new shortwave to store 
     186      IF( llnew ) THEN 
    170187          IF( lwp ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, & 
    171              &                      ' time = ', (iseclast+rdt*nn_dttrc/2.)/3600.,'hours ' 
    172           isecfst = iseclast 
     188             &                      ' time = ', rseclast/3600.,'hours ' 
     189          rsecfst = rseclast 
    173190          DO jn = 1, nb_rec_per_day - 1 
    174191             qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1) 
     
    182199         IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file  kt =', kt 
    183200         IF(lwp) WRITE(numout,*) '~~~~~~~' 
     201         zkt = REAL( kt, wp ) 
     202         CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt ) 
     203          DO jn = 1, nb_rec_per_day  
     204             IF( jn <= 9 )  THEN 
     205               WRITE(cl1,'(i1)') jn 
     206               CALL iom_rstput( kt, nitrst, numrtw, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) ) 
     207             ELSE 
     208               WRITE(cl2,'(i2.2)') jn 
     209               CALL iom_rstput( kt, nitrst, numrtw, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) ) 
     210             ENDIF 
     211         ENDDO 
    184212         CALL iom_rstput( kt, nitrst, numrtw, 'qsr_mean', qsr_mean(:,:) ) 
    185213      ENDIF 
    186      ! 
     214      ! 
    187215   END SUBROUTINE trc_mean_qsr 
    188216 
Note: See TracChangeset for help on using the changeset viewer.