Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90	(revision 5602)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90	(revision 6101)
@@ -27,4 +27,10 @@
    USE iom
    USE zdfmxl_crs
+   USE eosbn2
+   USE zdfevd_crs
+   USE zdftke
+   USE zdftke_crs
+
+!   USE ieee_arithmetic
 
    IMPLICIT NONE
@@ -40,5 +46,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
+   !! $Id $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -65,12 +71,18 @@
       !!
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zfse3t, zfse3u, zfse3v, zfse3w ! 3D workspace for e3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt, zs 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt, zs , ztmp
       REAL(wp), POINTER, DIMENSION(:,:)   :: z2d,z2d_crs
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_crs, zs_crs !
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_crs, zs_crs, zerr_crs,zmax_crs
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmp_crs
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: avte_crs
       REAL(wp)       :: z2dcrsu, z2dcrsv
-      REAL(wp)       :: zmin,zmax
+      REAL(wp)       :: zmin,zmax,icnt1,icnt2
       INTEGER :: i,j,ijis,ijie,ijjs,ijje
       REAL(wp)       :: zw,zwp1,zum1,zu,zvm1,zv,zsnm,zsm,z
+      REAL(wp)       :: zerr, zerr0, zerr1, zmean
+      INTEGER,DIMENSION(4,3) :: ind
+      REAL(wp),DIMENSION(4) :: zwgt
       INTEGER ::  iji,ijj
+      INTEGER :: jl,jm,jn
       !!
       !!----------------------------------------------------------------------
@@ -81,8 +93,10 @@
       CALL wrk_alloc( jpi, jpj, jpk, zfse3t, zfse3w )
       CALL wrk_alloc( jpi, jpj, jpk, zfse3u, zfse3v )
-      CALL wrk_alloc( jpi, jpj, jpk, zt, zs         )
+      CALL wrk_alloc( jpi, jpj, jpk, zt, zs , ztmp        )
       CALL wrk_alloc( jpi, jpj,      z2d            )
       !
-      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs )
+      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs, zerr_crs,zmax_crs )
+      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, ztmp_crs )
+      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, 4, avte_crs )
       CALL wrk_alloc( jpi_crs, jpj_crs, z2d_crs     )
 
@@ -129,5 +143,10 @@
       CALL iom_put( "sst" , tsn_crs(:,:,1,jp_tem) )    ! sst
 
-      
+      !n2 before
+      zt(:,:,:) = rn2b(:,:,:)  ;      zt_crs(:,:,:) = 0._wp
+      CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
+      rb2_crs(:,:,:) = zt_crs(:,:,:)
+      CALL iom_put("rb2_crs",rb2_crs)
+
       !  Salinity
       zs(:,:,:) = tsb(:,:,:,jp_sal)  ;      zs_crs(:,:,:) = 0._wp
@@ -252,4 +271,148 @@
          CASE ( 2 )
             CALL crs_dom_ope( avt, 'MIN', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
+         CASE ( 3 )
+            CALL crs_dom_ope( avt, 'LOGVOL', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, p_mask_crs=tmask_crs, psgn=1.0 )
+         CASE ( 4 )
+            CALL crs_dom_ope( avt, 'MED', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
+         CASE ( 5 )
+            CALL crs_dom_ope( en , 'VOL', 'W', tmask, en_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
+            CALL crs_dom_ope( taum , 'SUM', 'T', tmask, taum_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
+            CALL crs_dom_ope( rn2(:,:,:), 'VOL', 'W', tmask, rn2_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
+            IF( kt==nit000 )CALL tke_avn_ini_crs
+            CALL tke_avn_crs
+            CALL zdf_evd_crs(kt)
+         CASE ( 6 )
+
+            avte_crs(:,:,:,:) = 0._wp
+            ztmp(:,:,:)=1.
+
+            zt(:,:,:) = 0._wp
+            zs(:,:,:) = 0._wp
+            DO jk=2,jpk 
+               WHERE( fse3w(:,:,jk) .NE. 0._wp) zs(:,:,jk) = ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) / fse3w(:,:,jk)
+               zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
+            ENDDO
+            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
+            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
+            zt_crs=tmask_crs*zt_crs
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
+            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,1) = zt_crs / zs_crs
+            zmin=MINVAL(avte_crs(:,:,:,1));zmax=MAXVAL(avte_crs(:,:,:,1));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
+
+            zt(:,:,:) = 0._wp
+            zs(:,:,:) = 0._wp
+            DO jk=2,jpk 
+               WHERE( fse3w(:,:,jk) .NE. 0._wp) zs(:,:,jk) = ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) / fse3w(:,:,jk)
+               zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
+            ENDDO
+            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax 
+            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax 
+            zt_crs=tmask_crs*zt_crs
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax
+            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,2) = zt_crs / zs_crs
+            zmin=MINVAL(avte_crs(:,:,:,2));zmax=MAXVAL(avte_crs(:,:,:,2));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax
+
+            zt(:,:,:) = 0._wp
+            zs(:,:,:) = 0._wp
+            DO jk=2,jpk
+                WHERE( fse3w(:,:,jk) .NE. 0._wp ) zs(:,:,jk)= ( rn_a0 * ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) +  &
+                                                                  &   rn_b0 * ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) )/ fse3w(:,:,jk) 
+                zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)                                                           
+            ENDDO
+            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax 
+            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax 
+            zt_crs=tmask_crs*zt_crs
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax
+            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,3) = zt_crs / zs_crs
+            zmin=MINVAL(avte_crs(:,:,:,3));zmax=MAXVAL(avte_crs(:,:,:,3));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax
+
+            zt(:,:,:) = 0._wp
+            zs(:,:,:) = 0._wp
+            DO jk=2,jpk
+                WHERE( fse3w(:,:,jk) .NE. 0._wp ) zs(:,:,jk)= ( rn_a0 * ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) -  &
+                                                                  &   rn_b0 * ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) )/ fse3w(:,:,jk) 
+                zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
+            ENDDO
+            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax 
+            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax 
+            zt_crs=tmask_crs*zt_crs
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax
+            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,4) = zt_crs / zs_crs
+            zmin=MINVAL(avte_crs(:,:,:,4));zmax=MAXVAL(avte_crs(:,:,:,4));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax
+
+            CALL iom_put( "avte_crs1", avte_crs(:,:,:,1) )   !  Kz
+            CALL iom_put( "avte_crs2", avte_crs(:,:,:,2) )   !  Kz
+            CALL iom_put( "avte_crs3", avte_crs(:,:,:,3) )   !  Kz
+            CALL iom_put( "avte_crs4", avte_crs(:,:,:,4) )   !  Kz
+!--------------------- 
+            CALL crs_dom_ope( avt, 'MED', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=zfse3w, p_mask_crs=tmask_crs, psgn=1.0 )
+!?            zmin=MINVAL(zs_crs*tmask_crs);zmax=MAXVAL(zs_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"logvol zs_crs*tmask ",zmin,zmax ; call flush(numout)
+            CALL iom_put( "zs_crs", zs_crs )   !  Kzlogvol
+!--------------------- ok
+
+            CALL crs_dom_ope( avt, 'VOL', 'W', tmask, zmax_crs, p_e12=e1e2t, p_e3=zfse3w,  psgn=1.0 )
+            WRITE(narea+200,*)"zmax_crs ",SHAPE(zmax_crs) ; call flush(narea+200)
+            CALL iom_put( "zmax_crs", zmax_crs )   !  Kzlogvol
+            zmin=MINVAL(zmax_crs*tmask_crs);zmax=MAXVAL(zmax_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"vol zmax_crs*tmask ",zmin,zmax ; call flush(numout)
+!-------------------------nok
+            avt_crs=zs_crs
+
+
+            zerr0=0.01
+
+            icnt1=0
+            icnt2=0
+
+            zt_crs(:,:,:)=0._wp
+            zerr_crs(:,:,:)=0._wp
+            DO ji=1,jpi_crs 
+            DO jj=1,jpj_crs 
+            DO jk=1,jpk
+
+ 
+!--------------
+               zwgt(1:4)=0._wp
+               DO jm=1,4 ; IF( avte_crs(ji,jj,jk,jm)  .GE. 0._wp .AND.  avte_crs(ji,jj,jk,jm)  .LE. zmax_crs(ji,jj,jk)  ) zwgt(jm) = 1._wp ; ENDDO
+!--------------
+               IF( SUM(zwgt(1:4)) .NE. 0._wp ) THEN  
+                  zmean = SUM( zwgt(1:4)*avte_crs(ji,jj,jk,1:4)) / SUM(zwgt(1:4) )
+                  zerr  = SQRT(SUM( zwgt(1:4)*(avte_crs(ji,jj,jk,1:4)-zmean)*(avte_crs(ji,jj,jk,1:4)-zmean) ) / SUM(zwgt(1:4) ) )
+               ELSE
+                  zmean=0._wp
+                 zerr=1.e10
+               ENDIF
+!--------------
+
+               zerr_crs(ji,jj,jk)=zerr
+
+               IF( zerr .LE. zerr0 .AND. zmean .GT. 0._wp )zt_crs(ji,jj,jk)=zmean
+               IF( zerr .LE. zerr0 .AND. zmean .GT. 0._wp )avt_crs(ji,jj,jk)=zmean
+
+               IF( tmask_crs(ji,jj,jk) == 1 ) icnt1=icnt1+1
+               IF( tmask_crs(ji,jj,jk) == 1 .AND.  zerr .LE. zerr0 .AND. zmean .GT. 0._wp ) icnt2=icnt2+1
+
+               !IF( ieee_is_nan(  zt_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANMEANEFF ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
+               !IF( ieee_is_nan(  zs_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANLOG ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
+               !IF( ieee_is_nan( avt_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANAVT ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
+            ENDDO
+            ENDDO
+            ENDDO
+            zmin=MINVAL(avt_crs);zmax=MAXVAL(avt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avt_crs ",zmin,zmax  ; call flush(numout)
+            zmin=MINVAL(avt_crs*tmask_crs);zmax=MAXVAL(avt_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avt_crs*tmask ",zmin,zmax  ; call flush(numout)
+
+            CALL mpp_sum(icnt1)
+            CALL mpp_sum(icnt2)
+            IF(lwp)WRITE(numout,*)"TOTO",kt,icnt1,icnt2
+            CALL iom_put( "zt_crs", zt_crs )   !  Kz
+            CALL iom_put( "zerr_crs", zerr_crs )   !  Kz
+
       END SELECT
       !
@@ -293,7 +456,9 @@
       CALL wrk_dealloc( jpi, jpj, jpk, zfse3t, zfse3w )
       CALL wrk_dealloc( jpi, jpj, jpk, zfse3u, zfse3v )
-      CALL wrk_dealloc( jpi, jpj, jpk, zt, zs         )
+      CALL wrk_dealloc( jpi, jpj, jpk, zt, zs, ztmp   )
       CALL wrk_dealloc( jpi, jpj, z2d                 )
-      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs )
+      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs, zerr_crs,zmax_crs )
+      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, ztmp_crs )
+      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, 4, avte_crs )
       CALL wrk_dealloc( jpi_crs, jpj_crs, z2d_crs     )
       !
