Changeset 6455 for branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC

Timestamp:

2016-04-08T10:57:55+02:00 (8 years ago)

Author:

aumont

Message:

upgrade to last revision of 3_6_stable

Location:

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC

Files:

• DIA/diawri.F90 (modified) (3 diffs)
• DOM/domvvl.F90 (modified) (1 diff)
• IOM/iom.F90 (modified) (10 diffs)
• LBC/mppini.F90 (modified) (5 diffs)
• LBC/mppini_2.h90 (modified) (6 diffs)
• LDF/ldfslp.F90 (modified) (1 diff)
• LDF/ldftra_oce.F90 (modified) (3 diffs)
• LDF/ldftra_substitute.h90 (modified) (4 diffs)
• SBC/albedo.F90 (modified) (7 diffs)
• SBC/sbc_ice.F90 (modified) (2 diffs)
• SBC/sbcblk_clio.F90 (modified) (1 diff)
• SBC/sbcblk_core.F90 (modified) (2 diffs)
• SBC/sbccpl.F90 (modified) (10 diffs)
• SBC/sbcice_lim.F90 (modified) (9 diffs)
• SBC/sbcmod.F90 (modified) (1 diff)
• TRA/traldf.F90 (modified) (2 diffs)
• TRA/traqsr.F90 (modified) (12 diffs)
• ZDF/zdfddm.F90 (modified) (1 diff)
• ZDF/zdfmxl.F90 (modified) (3 diffs)
• ZDF/zdftke.F90 (modified) (2 diffs)
• ZDF/zdftmx.F90 (modified) (1 diff)
• step.F90 (modified) (1 diff)

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

@@ -145 +145 @@
       ENDIF
 
-      IF( .NOT.lk_vvl ) THEN
-         CALL iom_put( "e3t" , fse3t_n(:,:,:) )
-         CALL iom_put( "e3u" , fse3u_n(:,:,:) )
-         CALL iom_put( "e3v" , fse3v_n(:,:,:) )
-         CALL iom_put( "e3w" , fse3w_n(:,:,:) )
-      ENDIF
+      ! Output of initial vertical scale factor
+      CALL iom_put("e3t_0", e3t_0(:,:,:) )
+      CALL iom_put("e3u_0", e3t_0(:,:,:) )
+      CALL iom_put("e3v_0", e3t_0(:,:,:) )
+      !
+      CALL iom_put( "e3t" , fse3t_n(:,:,:) )
+      CALL iom_put( "e3u" , fse3u_n(:,:,:) )
+      CALL iom_put( "e3v" , fse3v_n(:,:,:) )
+      CALL iom_put( "e3w" , fse3w_n(:,:,:) )
+      IF( iom_use("e3tdef") )   &
+         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
+
 
       CALL iom_put( "ssh" , sshn )                 ! sea surface height
-      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
       
       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
@@ -243 +248 @@
       CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
       CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
+                                                            ! Log of eddy diff coef
+      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt  (:,:,:) ) ) )
+      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) )
 
       IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
@@ -307 +315 @@
          CALL iom_put( "eken", rke )           
       ENDIF
-         
+      !
+      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
+      !
       IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
          z3d(:,:,jpk) = 0.e0

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

@@ -665 +665 @@
          ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
       END DO
-
-      ! Write outputs
-      ! =============
-      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) )
-      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) )
-      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) )
-      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) )
-      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) )
-      IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
 
       ! write restart file

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

@@ -139 +139 @@
       ! horizontal grid definition
 
-#if ! defined key_xios2
       CALL set_scalar
-#endif
 
       IF( TRIM(cdname) == TRIM(cxios_context) ) THEN  
@@ -1193 +1191 @@
       REAL(wp), DIMENSION(:)   , OPTIONAL, INTENT(in) ::   lonvalue, latvalue
       REAL(wp), DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   bounds_lon, bounds_lat, area
-      LOGICAL,  DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   mask
+#if ! defined key_xios2
+     LOGICAL,  DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   mask
+#else
+      LOGICAL,  DIMENSION(:) , OPTIONAL, INTENT(in) ::   mask
+#endif
 
 #if ! defined key_xios2
@@ -1215 +1217 @@
          CALL xios_set_domain_attr     ( cdid, ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj,   &
             &    data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj ,   &
-            &    lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_2D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon,      &
+            &    lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_1D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon,                  &
             &    bounds_lat_1D=bounds_lat, area=area, type='curvilinear')
      ENDIF
@@ -1221 +1223 @@
          CALL xios_set_domaingroup_attr( cdid, ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj,   &
             &    data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj ,   &
-            &    lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_2D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon,      &
+            &    lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_1D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon,                  &
             &    bounds_lat_1D=bounds_lat, area=area, type='curvilinear' )
       ENDIF
@@ -1234 +1236 @@
      INTEGER                  , OPTIONAL, INTENT(in) ::   ibegin, jbegin, ni, nj
 
-     IF ( xios_is_valid_domain     (cdid) ) THEN
+     IF ( xios_is_valid_zoom_domain     (cdid) ) THEN
          CALL xios_set_zoom_domain_attr     ( cdid, ibegin=ibegin, jbegin=jbegin, ni=ni,    &
            &   nj=nj)
@@ -1326 +1328 @@
       IF ( xios_is_valid_gridgroup(cdid) )   CALL xios_set_gridgroup_attr( cdid, mask=mask )
 #else
-      IF ( xios_is_valid_grid     (cdid) )   CALL xios_set_grid_attr     ( cdid, mask3=mask )
-      IF ( xios_is_valid_gridgroup(cdid) )   CALL xios_set_gridgroup_attr( cdid, mask3=mask )
+      IF ( xios_is_valid_grid     (cdid) )   CALL xios_set_grid_attr     ( cdid, mask_3D=mask )
+      IF ( xios_is_valid_gridgroup(cdid) )   CALL xios_set_gridgroup_attr( cdid, mask_3D=mask )
 #endif
       CALL xios_solve_inheritance()
@@ -1388 +1390 @@
          END SELECT
          !
+#if ! defined key_xios2
          CALL iom_set_domain_attr( "grid_"//cdgrd       , mask = RESHAPE(zmask(nldi:nlei,nldj:nlej,1),(/ni,nj    /)) /= 0. )
+#else
+         CALL iom_set_domain_attr( "grid_"//cdgrd       , mask = RESHAPE(zmask(nldi:nlei,nldj:nlej,1),(/ni*nj    /)) /= 0. )
+#endif  
          CALL iom_set_grid_attr  ( "grid_"//cdgrd//"_3D", mask = RESHAPE(zmask(nldi:nlei,nldj:nlej,:),(/ni,nj,jpk/)) /= 0. )
       ENDIF
@@ -1532 +1538 @@
 #else
 ! Pas teste : attention aux indices !
-      CALL iom_set_domain_attr("ptr", ni_glo=jpiglo, nj_glo=jpjglo, ibegin=nimpp+nldi-2, jbegin=njmpp+nldj-2, ni=ni, nj=nj)
-      CALL iom_set_domain_attr("ptr", data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj)
-      CALL iom_set_domain_attr("ptr", lonvalue = zlon,   &
+      CALL iom_set_domain_attr("gznl", ni_glo=jpiglo, nj_glo=jpjglo, ibegin=nimpp+nldi-2, jbegin=njmpp+nldj-2, ni=ni, nj=nj)
+      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj)
+      CALL iom_set_domain_attr("gznl", lonvalue = zlon,   &
          &                             latvalue = RESHAPE(plat(nldi:nlei, nldj:nlej),(/ ni*nj /)))  
-       CALL iom_set_zoom_domain_attr ('ptr', ibegin=ix, nj=jpjglo)
+       CALL iom_set_zoom_domain_attr ("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=jpjglo)
 #endif
 
@@ -1552 +1558 @@
       REAL(wp), DIMENSION(1)   ::   zz = 1.
       !!----------------------------------------------------------------------
+#if ! defined key_xios2
       CALL iom_set_domain_attr('scalarpoint', ni_glo=jpnij, nj_glo=1, ibegin=narea, jbegin=1, ni=1, nj=1)
+#else
+      CALL iom_set_domain_attr('scalarpoint', ni_glo=jpnij, nj_glo=1, ibegin=narea-1, jbegin=0, ni=1, nj=1)
+#endif
       CALL iom_set_domain_attr('scalarpoint', data_dim=2, data_ibegin = 1, data_ni = 1, data_jbegin = 1, data_nj = 1)
       
@@ -1778 +1788 @@
             idx = INDEX(clname,'@freq@') + INDEX(clname,'@FREQ@')
             DO WHILE ( idx /= 0 ) 
-              IF ( output_freq%hour /= 0 ) THEN
+              IF ( output_freq%timestep /= 0) THEN
+                  WRITE(clfreq,'(I18,A2)')INT(output_freq%timestep),'ts' 
+                  itrlen = LEN_TRIM(ADJUSTL(clfreq))
+              ELSE IF ( output_freq%hour /= 0 ) THEN
                   WRITE(clfreq,'(I19,A1)')INT(output_freq%hour),'h' 
                   itrlen = LEN_TRIM(ADJUSTL(clfreq))

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90

@@ -201 +201 @@
       
 #endif
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '           defines mpp subdomains'
-         WRITE(numout,*) '           ----------------------'
-         WRITE(numout,*) '           iresti=',iresti,' irestj=',irestj
-         WRITE(numout,*) '           jpni  =',jpni  ,' jpnj  =',jpnj
-         ifreq = 4
-         il1   = 1
-         DO jn = 1, (jpni-1)/ifreq+1
-            il2 = MIN( jpni, il1+ifreq-1 )
-            WRITE(numout,*)
-            WRITE(numout,9200) ('***',ji = il1,il2-1)
-            DO jj = jpnj, 1, -1
-               WRITE(numout,9203) ('   ',ji = il1,il2-1)
-               WRITE(numout,9202) jj, ( ilcit(ji,jj),ilcjt(ji,jj),ji = il1,il2 )
-               WRITE(numout,9203) ('   ',ji = il1,il2-1)
-               WRITE(numout,9200) ('***',ji = il1,il2-1)
-            END DO
-            WRITE(numout,9201) (ji,ji = il1,il2)
-            il1 = il1+ifreq
-         END DO
- 9200    FORMAT('     ***',20('*************',a3))
- 9203    FORMAT('     *     ',20('         *   ',a3))
- 9201    FORMAT('        ',20('   ',i3,'          '))
- 9202    FORMAT(' ',i3,' *  ',20(i3,'  x',i3,'   *   '))
-      ENDIF
-
-      zidom = nreci
-      DO ji = 1, jpni
-         zidom = zidom + ilcit(ji,1) - nreci
-      END DO
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*)' sum ilcit(i,1) = ', zidom, ' jpiglo = ', jpiglo
-      
-      zjdom = nrecj
-      DO jj = 1, jpnj
-         zjdom = zjdom + ilcjt(1,jj) - nrecj
-      END DO
-      IF(lwp) WRITE(numout,*)' sum ilcit(1,j) = ', zjdom, ' jpjglo = ', jpjglo
-      IF(lwp) WRITE(numout,*)
-      
 
       !  2. Index arrays for subdomains
@@ -304 +263 @@
          nlejt(jn) = nlej
       END DO
-      
-
-      ! 4. From global to local
+
+      ! 4. Subdomain print
+      ! ------------------
+      
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' mpp_init: defines mpp subdomains'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~  ----------------------'
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'iresti=',iresti,' irestj=',irestj
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'jpni=',jpni,' jpnj=',jpnj
+      zidom = nreci
+      DO ji = 1, jpni
+         zidom = zidom + ilcit(ji,1) - nreci
+      END DO
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)' sum ilcit(i,1)=', zidom, ' jpiglo=', jpiglo
+
+      zjdom = nrecj
+      DO jj = 1, jpnj
+         zjdom = zjdom + ilcjt(1,jj) - nrecj
+      END DO
+      IF(lwp) WRITE(numout,*)' sum ilcit(1,j)=', zjdom, ' jpjglo=', jpjglo
+      IF(lwp) WRITE(numout,*)
+
+      IF(lwp) THEN
+         ifreq = 4
+         il1   = 1
+         DO jn = 1, (jpni-1)/ifreq+1
+            il2 = MIN( jpni, il1+ifreq-1 )
+            WRITE(numout,*)
+            WRITE(numout,9200) ('***',ji = il1,il2-1)
+            DO jj = jpnj, 1, -1
+               WRITE(numout,9203) ('   ',ji = il1,il2-1)
+               WRITE(numout,9202) jj, ( ilcit(ji,jj),ilcjt(ji,jj),ji = il1,il2 )
+               WRITE(numout,9204) (nfipproc(ji,jj),ji=il1,il2)
+               WRITE(numout,9203) ('   ',ji = il1,il2-1)
+               WRITE(numout,9200) ('***',ji = il1,il2-1)
+            END DO
+            WRITE(numout,9201) (ji,ji = il1,il2)
+            il1 = il1+ifreq
+         END DO
+ 9200     FORMAT('     ***',20('*************',a3))
+ 9203     FORMAT('     *     ',20('         *   ',a3))
+ 9201     FORMAT('        ',20('   ',i3,'          '))
+ 9202     FORMAT(' ',i3,' *  ',20(i3,'  x',i3,'   *   '))
+ 9204     FORMAT('     *  ',20('      ',i3,'   *   '))
+      ENDIF
+
+      ! 5. From global to local
       ! -----------------------
 
@@ -313 +319 @@
 
 
-      ! 5. Subdomain neighbours
+      ! 6. Subdomain neighbours
       ! ----------------------
 
@@ -436 +442 @@
          WRITE(numout,*) ' nimpp  = ', nimpp
          WRITE(numout,*) ' njmpp  = ', njmpp
-         WRITE(numout,*) ' nbse   = ', nbse  , ' npse   = ', npse
-         WRITE(numout,*) ' nbsw   = ', nbsw  , ' npsw   = ', npsw
-         WRITE(numout,*) ' nbne   = ', nbne  , ' npne   = ', npne
-         WRITE(numout,*) ' nbnw   = ', nbnw  , ' npnw   = ', npnw
+         WRITE(numout,*) ' nreci  = ', nreci  , ' npse   = ', npse
+         WRITE(numout,*) ' nrecj  = ', nrecj  , ' npsw   = ', npsw
+         WRITE(numout,*) ' jpreci = ', jpreci , ' npne   = ', npne
+         WRITE(numout,*) ' jprecj = ', jprecj , ' npnw   = ', npnw
+         WRITE(numout,*)
       ENDIF
 
@@ -446 +453 @@
       ! Prepare mpp north fold
 
-      IF (jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
+      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
          CALL mpp_ini_north
-      END IF
+         IF(lwp) WRITE(numout,*) ' mpp_init : North fold boundary prepared for jpni >1'
+      ENDIF
 
       ! Prepare NetCDF output file (if necessary)

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90

@@ -318 +318 @@
          ENDIF
 
+         ! Check wet points over the entire domain to preserve the MPI communication stencil
          isurf = 0
-         DO jj = 1+jprecj, ilj-jprecj
-            DO  ji = 1+jpreci, ili-jpreci
+         DO jj = 1, ilj
+            DO  ji = 1, ili
                IF( imask(ji+iimppt(ii,ij)-1, jj+ijmppt(ii,ij)-1) == 1) isurf = isurf+1
             END DO
          END DO
+
          IF(isurf /= 0) THEN
             icont = icont + 1
@@ -333 +335 @@
 
       nfipproc(:,:) = ipproc(:,:)
-
 
       ! Control
@@ -441 +442 @@
       ii = iin(narea)
       ij = ijn(narea)
+
+      ! set default neighbours
+      noso = ioso(ii,ij)
+      nowe = iowe(ii,ij)
+      noea = ioea(ii,ij)
+      nono = iono(ii,ij) 
+      npse = iose(ii,ij)
+      npsw = iosw(ii,ij)
+      npne = ione(ii,ij)
+      npnw = ionw(ii,ij)
+
+      ! check neighbours location
       IF( ioso(ii,ij) >= 0 .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN 
          iiso = 1 + MOD(ioso(ii,ij),jpni)
@@ -511 +524 @@
       IF (lwp) THEN
          CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
+         WRITE(inum,'(a)') '   jpnij     jpi     jpj     jpk  jpiglo  jpjglo'
          WRITE(inum,'(6i8)') jpnij,jpi,jpj,jpk,jpiglo,jpjglo
          WRITE(inum,'(a)') 'NAREA nlci nlcj nldi nldj nlei nlej nimpp njmpp'
@@ -523 +537 @@
       END IF
 
-      IF( nperio == 1 .AND.jpni /= 1 ) CALL ctl_stop( ' mpp_init2:  error on cyclicity' )
-
-      ! Prepare mpp north fold
-
-      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
-         CALL mpp_ini_north
-         IF(lwp) WRITE(numout,*) ' mpp_init2 : North fold boundary prepared for jpni >1'
-      ENDIF
-
       ! Defined npolj, either 0, 3 , 4 , 5 , 6
       ! In this case the important thing is that npolj /= 0
@@ -548 +553 @@
       ENDIF
 
+      ! Periodicity : no corner if nbondi = 2 and nperio != 1
+
+      IF(lwp) THEN
+         WRITE(numout,*) ' nproc  = ', nproc
+         WRITE(numout,*) ' nowe   = ', nowe  , ' noea   =  ', noea
+         WRITE(numout,*) ' nono   = ', nono  , ' noso   =  ', noso
+         WRITE(numout,*) ' nbondi = ', nbondi
+         WRITE(numout,*) ' nbondj = ', nbondj
+         WRITE(numout,*) ' npolj  = ', npolj
+         WRITE(numout,*) ' nperio = ', nperio
+         WRITE(numout,*) ' nlci   = ', nlci
+         WRITE(numout,*) ' nlcj   = ', nlcj
+         WRITE(numout,*) ' nimpp  = ', nimpp
+         WRITE(numout,*) ' njmpp  = ', njmpp
+         WRITE(numout,*) ' nreci  = ', nreci  , ' npse   = ', npse
+         WRITE(numout,*) ' nrecj  = ', nrecj  , ' npsw   = ', npsw
+         WRITE(numout,*) ' jpreci = ', jpreci , ' npne   = ', npne
+         WRITE(numout,*) ' jprecj = ', jprecj , ' npnw   = ', npnw
+         WRITE(numout,*)
+      ENDIF
+
+      IF( nperio == 1 .AND. jpni /= 1 ) CALL ctl_stop( ' mpp_init2: error on cyclicity' )
+
+      ! Prepare mpp north fold
+
+      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
+         CALL mpp_ini_north
+         IF(lwp) WRITE(numout,*) ' mpp_init2 : North fold boundary prepared for jpni >1'
+      ENDIF
+
       ! Prepare NetCDF output file (if necessary)
       CALL mpp_init_ioipsl
 
-      ! Periodicity : no corner if nbondi = 2 and nperio != 1
-
-      IF(lwp) THEN
-         WRITE(numout,*) ' nproc=  ',nproc
-         WRITE(numout,*) ' nowe=   ',nowe
-         WRITE(numout,*) ' noea=   ',noea
-         WRITE(numout,*) ' nono=   ',nono
-         WRITE(numout,*) ' noso=   ',noso
-         WRITE(numout,*) ' nbondi= ',nbondi
-         WRITE(numout,*) ' nbondj= ',nbondj
-         WRITE(numout,*) ' npolj=  ',npolj
-         WRITE(numout,*) ' nperio= ',nperio
-         WRITE(numout,*) ' nlci=   ',nlci
-         WRITE(numout,*) ' nlcj=   ',nlcj
-         WRITE(numout,*) ' nimpp=  ',nimpp
-         WRITE(numout,*) ' njmpp=  ',njmpp
-         WRITE(numout,*) ' nbse=   ',nbse,' npse= ',npse
-         WRITE(numout,*) ' nbsw=   ',nbsw,' npsw= ',npsw
-         WRITE(numout,*) ' nbne=   ',nbne,' npne= ',npne
-         WRITE(numout,*) ' nbnw=   ',nbnw,' npnw= ',npnw
-      ENDIF
 
    END SUBROUTINE mpp_init2

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

@@ -188 +188 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
-                  IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp)
-                  IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp)
-                  IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp)
-                  IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp)
-                  IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp)
+               zhmlpu(ji,jj) = ( MAX(hmlpt(ji,jj)  , hmlpt  (ji+1,jj  ), 5._wp)   &
+                  &            - MAX(risfdep(ji,jj), risfdep(ji+1,jj  )       )   )
+               zhmlpv(ji,jj) = ( MAX(hmlpt  (ji,jj), hmlpt  (ji  ,jj+1), 5._wp)   &
+                  &            - MAX(risfdep(ji,jj), risfdep(ji  ,jj+1)       )   )
                ENDDO
             ENDDO

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90

@@ -41 +41 @@
 
    REAL(wp), PUBLIC ::   rldf                        !: multiplicative factor of diffusive coefficient
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   r_fact_lap
                                                      !: Needed to define the ratio between passive and active tracer diffusion coef. 
 
@@ -92 +93 @@
       !!                 ***  FUNCTION ldftra_oce_alloc  ***
      !!----------------------------------------------------------------------
-     INTEGER, DIMENSION(3) :: ierr
+     INTEGER, DIMENSION(4) :: ierr
      !!----------------------------------------------------------------------
      ierr(:) = 0
@@ -116 +117 @@
 # endif
 #endif
+      ALLOCATE( r_fact_lap(jpi,jpj,jpk), STAT=ierr(4) )
       ldftra_oce_alloc = MAXVAL( ierr )
       IF( ldftra_oce_alloc /= 0 )   CALL ctl_warn('ldftra_oce_alloc: failed to allocate arrays')

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_substitute.h90

@@ -13 +13 @@
 !   'key_traldf_c3d' :                 aht: 3D coefficient
 #       define   fsahtt(i,j,k)   rldf * ahtt(i,j,k)
-#       define   fsahtu(i,j,k)   rldf * ahtu(i,j,k)
+#       define   fsahtu(i,j,k)   rldf * ahtu(i,j,k) * r_fact_lap(i,j,k)
 #       define   fsahtv(i,j,k)   rldf * ahtv(i,j,k)
 #       define   fsahtw(i,j,k)   rldf * ahtw(i,j,k)
@@ -19 +19 @@
 !   'key_traldf_c2d' :                 aht: 2D coefficient
 #       define   fsahtt(i,j,k)   rldf * ahtt(i,j)
-#       define   fsahtu(i,j,k)   rldf * ahtu(i,j)
+#       define   fsahtu(i,j,k)   rldf * ahtu(i,j) * r_fact_lap(i,j,k)
 #       define   fsahtv(i,j,k)   rldf * ahtv(i,j)
 #       define   fsahtw(i,j,k)   rldf * ahtw(i,j)
@@ -25 +25 @@
 !   'key_traldf_c1d' :                aht: 1D coefficient
 #       define   fsahtt(i,j,k)   rldf * ahtt(k)
-#       define   fsahtu(i,j,k)   rldf * ahtu(k)
+#       define   fsahtu(i,j,k)   rldf * ahtu(k) * r_fact_lap(i,j,k)
 #       define   fsahtv(i,j,k)   rldf * ahtv(k)
 #       define   fsahtw(i,j,k)   rldf * ahtw(k)
@@ -31 +31 @@
 !   Default option :             aht: Constant coefficient
 #      define   fsahtt(i,j,k)   rldf * aht0
-#      define   fsahtu(i,j,k)   rldf * aht0
+#      define   fsahtu(i,j,k)   rldf * aht0 * r_fact_lap(i,j,k)
 #      define   fsahtv(i,j,k)   rldf * aht0
 #      define   fsahtw(i,j,k)   rldf * aht0

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

@@ -9 +9 @@
    !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0
    !!             -   ! 2006-08  (G. Madec)  cleaning for surface module
+   !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo
    !!----------------------------------------------------------------------
 
@@ -29 +30 @@
 
    INTEGER  ::   albd_init = 0      !: control flag for initialization
-   REAL(wp) ::   zzero     = 0.e0   ! constant values
-   REAL(wp) ::   zone      = 1.e0   !    "       "
-
-   REAL(wp) ::   c1     = 0.05    ! constants values
-   REAL(wp) ::   c2     = 0.10    !    "        "
-   REAL(wp) ::   rmue   = 0.40    !  cosine of local solar altitude
-
+  
+   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude
+   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
+   REAL(wp) ::   c1       = 0.05    ! snow thickness (only for nn_ice_alb=0)
+   REAL(wp) ::   c2       = 0.10    !  "        "
+   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
+ 
    !                             !!* namelist namsbc_alb
-   REAL(wp) ::   rn_cloud         !  cloudiness effect on snow or ice albedo (Grenfell & Perovich, 1984)
-#if defined key_lim3
-   REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-#else
-   REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-#endif
-   REAL(wp) ::   rn_alphd         !  coefficients for linear interpolation used to compute
-   REAL(wp) ::   rn_alphdi        !  albedo between two extremes values (Pyane, 1972)
-   REAL(wp) ::   rn_alphc         ! 
+   INTEGER  ::   nn_ice_alb
+   REAL(wp) ::   rn_albice
 
    !!----------------------------------------------------------------------
@@ -59 +53 @@
       !!          
       !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!                as well as the ocean one
       !!       
-      !! ** Method  : - Computation of the albedo of snow or ice (choose the 
-      !!                rignt one by a large number of tests
-      !!              - Computation of the albedo of the ocean
-      !!
-      !! References :   Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
+      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
+      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
+      !!                                                                           and Grenfell & Perovich (JGR 2004)
+      !!                Description of scheme 1:
+      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
+      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
+      !!                     0-5cm  : linear function of ice thickness
+      !!                     5-150cm: log    function of ice thickness
+      !!                     > 150cm: constant
+      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
+      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
+      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
+      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
+      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
+      !!
+      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
+      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
+      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
+      !!                     under melting conditions than under freezing conditions
+      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
+      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
+      !!
+      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !!                Brandt et al. 2005, J. Climate, vol 18
+      !!                Grenfell & Perovich 2004, JGR, vol 109 
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
@@ -73 +87 @@
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky
       !!
-      INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      INTEGER  ::   ijpl          ! number of ice categories (3rd dim of ice input arrays)
-      REAL(wp) ::   zalbpsnm      ! albedo of ice under clear sky when snow is melting
-      REAL(wp) ::   zalbpsnf      ! albedo of ice under clear sky when snow is freezing
-      REAL(wp) ::   zalbpsn       ! albedo of snow/ice system when ice is coverd by snow
-      REAL(wp) ::   zalbpic       ! albedo of snow/ice system when ice is free of snow
-      REAL(wp) ::   zithsn        ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow)
-      REAL(wp) ::   zitmlsn       ! = 1 freezinz snow (pt_ice >=rt0_snow) ; = 0 melting snow (pt_ice<rt0_snow)
-      REAL(wp) ::   zihsc1        ! = 1 hsn <= c1 ; = 0 hsn > c1
-      REAL(wp) ::   zihsc2        ! = 1 hsn >= c2 ; = 0 hsn < c2
-      !!
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalbfz    ! = rn_alphdi for freezing ice ; = rn_albice for melting ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zficeth   !  function of ice thickness
+      INTEGER  ::   ji, jj, jl         ! dummy loop indices
+      INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays)
+      REAL(wp)            ::   ralb_im, ralb_sf, ralb_sm, ralb_if
+      REAL(wp)            ::   zswitch, z1_c1, z1_c2
+      REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free)
       !!---------------------------------------------------------------------
-      
+
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
-
-      CALL wrk_alloc( jpi,jpj,ijpl, zalbfz, zficeth )
+      
+      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
 
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
 
-      !---------------------------
-      !  Computation of  zficeth
-      !---------------------------
-      ! ice free of snow and melts
-      WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice
-      ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi
-      END  WHERE
-
-      WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz
-      ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 )
-      ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              &
-         &                                                                 - 0.8608 * ph_ice * ph_ice     &
-         &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
-      ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice
-      END WHERE
-
-!!gm old code
-!      DO jl = 1, ijpl
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN
-!                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl)
-!               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN
-!                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 )
-!               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN
-!                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               &
-!                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 &
-!                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl)
-!               ELSE
-!                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl) 
-!               ENDIF
-!            END DO
-!         END DO
-!      END DO
-!!gm end old code
-      
-      !----------------------------------------------- 
-      !    Computation of the snow/ice albedo system 
-      !-------------------------- ---------------------
-      
-      !    Albedo of snow-ice for clear sky.
-      !-----------------------------------------------    
-      DO jl = 1, ijpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !  Case of ice covered by snow.             
-               !                                        !  freezing snow        
-               zihsc1   = 1.0 - MAX( zzero , SIGN( zone , - ( ph_snw(ji,jj,jl) - c1 ) ) )
-               zalbpsnf = ( 1.0 - zihsc1 ) * (  zficeth(ji,jj,jl)                                             &
-                  &                           + ph_snw(ji,jj,jl) * ( rn_alphd - zficeth(ji,jj,jl) ) / c1  )   &
-                  &     +         zihsc1   * rn_alphd  
-               !                                        !  melting snow                
-               zihsc2   = MAX( zzero , SIGN( zone , ph_snw(ji,jj,jl) - c2 ) )
-               zalbpsnm = ( 1.0 - zihsc2 ) * ( rn_albice + ph_snw(ji,jj,jl) * ( rn_alphc - rn_albice ) / c2 )   &
-                  &     +         zihsc2   *   rn_alphc 
-               !
-               zitmlsn  =  MAX( zzero , SIGN( zone , pt_ice(ji,jj,jl) - rt0_snow ) )   
-               zalbpsn  =  zitmlsn * zalbpsnm + ( 1.0 - zitmlsn ) * zalbpsnf
-            
-               !  Case of ice free of snow.
-               zalbpic  = zficeth(ji,jj,jl) 
-            
-               ! albedo of the system   
-               zithsn   = 1.0 - MAX( zzero , SIGN( zone , - ph_snw(ji,jj,jl) ) )
-               pa_ice_cs(ji,jj,jl) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic
+      
+      SELECT CASE ( nn_ice_alb )
+
+      !------------------------------------------
+      !  Shine and Henderson-Sellers (1985)
+      !------------------------------------------
+      CASE( 0 )
+       
+         ralb_sf = 0.80       ! dry snow
+         ralb_sm = 0.65       ! melting snow
+         ralb_if = 0.72       ! bare frozen ice
+         ralb_im = rn_albice  ! bare puddled ice 
+         
+         !  Computation of ice albedo (free of snow)
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
+         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
+         END  WHERE
+      
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
+            &                                                                 - 0.8608 * ph_ice * ph_ice     &
+            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
+         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
+         END WHERE
+     
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! freezing snow
+                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
+                  !                                        !  freezing snow        
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
+                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
+                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
+                     &        +         zswitch   * ralb_sf  
+
+                  ! melting snow
+                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
+                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
+                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
+                      &     +         zswitch   *   ralb_sm 
+                  !
+                  ! snow albedo
+                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+               
+                  ! Ice/snow albedo
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
+                  !
+               END DO
             END DO
          END DO
-      END DO
-      
-      !    Albedo of snow-ice for overcast sky.
-      !----------------------------------------------  
-      pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rn_cloud       ! Oberhuber correction
-      !
-      CALL wrk_dealloc( jpi,jpj,ijpl, zalbfz, zficeth )
+
+         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
+
+      !------------------------------------------
+      !  New parameterization (2016)
+      !------------------------------------------
+      CASE( 1 ) 
+
+         ralb_im = rn_albice  ! bare puddled ice
+! compilation of values from literature
+         ralb_sf = 0.85      ! dry snow
+         ralb_sm = 0.75      ! melting snow
+         ralb_if = 0.60      ! bare frozen ice
+! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
+!         ralb_sf = 0.85       ! dry snow
+!         ralb_sm = 0.72       ! melting snow
+!         ralb_if = 0.65       ! bare frozen ice
+! Brandt et al 2005 (East Antarctica)
+!         ralb_sf = 0.87      ! dry snow
+!         ralb_sm = 0.82      ! melting snow
+!         ralb_if = 0.54      ! bare frozen ice
+! 
+         !  Computation of ice albedo (free of snow)
+         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+         z1_c2 = 1. / 0.05
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
+         ELSE WHERE                                              ;   zalb = ralb_if
+         END  WHERE
+         
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
+            &                                                                     ( LOG(1.5) - LOG(ph_ice) )
+         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice
+         END WHERE
+
+         z1_c1 = 1. / 0.02
+         z1_c2 = 1. / 0.03
+         !  Computation of the snow/ice albedo
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
+                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
+
+                   ! snow albedo
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+
+                  ! Ice/snow albedo   
+                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
+
+              END DO
+            END DO
+         END DO
+         ! Effect of the clouds (2d order polynomial)
+         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
+
+      END SELECT
+      
+      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
       !
    END SUBROUTINE albedo_ice
@@ -181 +229 @@
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
       !!
-      REAL(wp) ::   zcoef   ! local scalar
-      !!----------------------------------------------------------------------
-      !
-      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )      ! Parameterization of Briegled and Ramanathan, 1982 
-      pa_oce_cs(:,:) = zcoef               
-      pa_oce_os(:,:)  = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972
+      REAL(wp) :: zcoef 
+      !!----------------------------------------------------------------------
+      !
+      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
+      pa_oce_cs(:,:) = zcoef 
+      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
       !
    END SUBROUTINE albedo_oce
@@ -200 +248 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ rn_cloud, rn_albice, rn_alphd, rn_alphdi, rn_alphc
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
       !!----------------------------------------------------------------------
       !
@@ -219 +267 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '   Namelist namsbc_alb : albedo '
-         WRITE(numout,*) '      correction for snow and ice albedo                  rn_cloud  = ', rn_cloud
-         WRITE(numout,*) '      albedo of melting ice in the arctic and antarctic   rn_albice = ', rn_albice
-         WRITE(numout,*) '      coefficients for linear                             rn_alphd  = ', rn_alphd
-         WRITE(numout,*) '      interpolation used to compute albedo                rn_alphdi = ', rn_alphdi
-         WRITE(numout,*) '      between two extremes values (Pyane, 1972)           rn_alphc  = ', rn_alphc
+         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
+         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice
       ENDIF
       !

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

@@ -80 +80 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_oce       !: heat flux of precip and evap over ocean     [W/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_ice       !: heat flux of precip and evap over ice       [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: heat flux of precip over ice                [J/m3]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qevap_ice      !: heat flux of evap over ice                  [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: enthalpy of precip over ice                 [J/m3]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s]
 #endif
@@ -144 +145 @@
 #endif
 #if defined key_lim3
-         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,  &
-         &      qemp_ice(jpi,jpj)     , qemp_oce(jpi,jpj)      ,                       &
-         &      qns_oce (jpi,jpj)     , qsr_oce (jpi,jpj)      , emp_oce (jpi,jpj)  ,  &
+         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   &
+         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   &
+         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce (jpi,jpj)  ,   &
 #endif
          &      emp_ice(jpi,jpj)      ,  STAT= ierr(1) )

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

@@ -684 +684 @@
       qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
 
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
+      DO jl = 1, jpl
+         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) - lfus )
+                                   ! but then qemp_ice should also include sublimation 
+      END DO
+
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
 #endif

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

@@ -612 +612 @@
       ! --- evaporation --- !
       z1_lsub = 1._wp / Lsub
-      evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation
-      devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub
-      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean
+      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation
+      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT
+      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean
 
       ! --- evaporation minus precipitation --- !
@@ -637 +637 @@
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
       qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
+      DO jl = 1, jpl
+         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
+                                   ! But we do not have Tice => consider it at 0°C => evap=0 
+      END DO
 
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

@@ -1378 +1378 @@
       !
       INTEGER ::   jl         ! dummy loop index
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk, zsnw
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
       !
-      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
-      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
+      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
+      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
+      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
 
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -1423 +1425 @@
       END SELECT
 
-      IF( iom_use('subl_ai_cea') )   &
-         CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
-      !   
-      !                                                           ! runoffs and calving (put in emp_tot)
+#if defined key_lim3
+      ! zsnw = snow percentage over ice after wind blowing
+      zsnw(:,:) = 0._wp
+      CALL lim_thd_snwblow( p_frld, zsnw )
+      
+      ! --- evaporation (kg/m2/s) --- !
+      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1)
+      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
+      ! therefore, sublimation is not redistributed over the ice categories in case no subgrid scale fluxes are provided by atm.
+      zdevap_ice(:,:) = 0._wp
+      
+      ! --- evaporation minus precipitation corrected for the effect of wind blowing on snow --- !
+      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) - zsprecip * (1._wp - zsnw)
+      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip * (1._wp - zsnw)          
+
+      ! Sublimation over sea-ice (cell average)
+      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:) )
+      ! runoffs and calving (put in emp_tot)
+      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+      IF( srcv(jpr_cal)%laction ) THEN 
+         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
+         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
+      ENDIF
+
+      IF( ln_mixcpl ) THEN
+         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
+         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
+         emp_oce(:,:) = emp_oce(:,:) * xcplmask(:,:,0) + zemp_oce(:,:) * zmsk(:,:)
+         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
+         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
+         DO jl=1,jpl
+            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:)
+            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:)
+         ENDDO
+      ELSE
+         emp_tot(:,:) =         zemp_tot(:,:)
+         emp_ice(:,:) =         zemp_ice(:,:)
+         emp_oce(:,:) =         zemp_oce(:,:)     
+         sprecip(:,:) =         zsprecip(:,:)
+         tprecip(:,:) =         ztprecip(:,:)
+         DO jl=1,jpl
+            evap_ice (:,:,jl) = zevap_ice (:,:)
+            devap_ice(:,:,jl) = zdevap_ice(:,:)
+         ENDDO
+      ENDIF
+
+                                     CALL iom_put( 'snowpre'    , sprecip                         )  ! Snow
+      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw ) )  ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw   )  ! Snow over sea-ice         (cell average)    
+#else
+      ! Sublimation over sea-ice (cell average)
+      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )
+      ! runoffs and calving (put in emp_tot)
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
       IF( srcv(jpr_cal)%laction ) THEN 
@@ -1450 +1501 @@
       IF( iom_use('snow_ai_cea') )   &
          CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
+#endif
 
       !                                                      ! ========================= !
@@ -1505 +1557 @@
       IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
 
-#if defined key_lim3
-      CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
-
+#if defined key_lim3      
       ! --- evaporation --- !
-      ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation
-      ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice
-      !                 but it is incoherent WITH the ice model  
-      DO jl=1,jpl
-         evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1)
-      ENDDO
       zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean
-
-      ! --- evaporation minus precipitation --- !
-      emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:)
 
       ! --- non solar flux over ocean --- !
@@ -1525 +1566 @@
       WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
 
-      ! --- heat flux associated with emp --- !
-      zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing
+      ! --- heat flux associated with emp (W/m2) --- !
       zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
          &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip
          &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean
-      qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
-         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
-
+!      zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
+!         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
+      zqemp_ice(:,:) =      zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice (only)
+                                                                                                       ! qevap_ice=0 since we consider Tice=0°C
+      
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
       zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
 
-      ! --- total non solar flux --- !
-      zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:)
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
+      DO jl = 1, jpl
+         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0°C
+      END DO
+
+      ! --- total non solar flux (including evap/precip) --- !
+      zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:)
 
       ! --- in case both coupled/forced are active, we must mix values --- ! 
@@ -1545 +1591 @@
          qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
          DO jl=1,jpl
-            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
+            qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:)
+            qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:)
          ENDDO
          qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
          qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
-!!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0)
+         qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:)
       ELSE
          qns_tot  (:,:  ) = zqns_tot  (:,:  )
          qns_oce  (:,:  ) = zqns_oce  (:,:  )
          qns_ice  (:,:,:) = zqns_ice  (:,:,:)
-         qprec_ice(:,:)   = zqprec_ice(:,:)
-         qemp_oce (:,:)   = zqemp_oce (:,:)
-      ENDIF
-
-      CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
+         qevap_ice(:,:,:) = zqevap_ice(:,:,:)
+         qprec_ice(:,:  ) = zqprec_ice(:,:  )
+         qemp_oce (:,:  ) = zqemp_oce (:,:  )
+         qemp_ice (:,:  ) = zqemp_ice (:,:  )
+      ENDIF
 #else
-
       ! clem: this formulation is certainly wrong... but better than it was...
       zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with:
@@ -1577 +1623 @@
          qns_ice(:,:,:) = zqns_ice(:,:,:)
       ENDIF
-
 #endif
 
@@ -1628 +1673 @@
 
 #if defined key_lim3
-      CALL wrk_alloc( jpi,jpj, zqsr_oce ) 
       ! --- solar flux over ocean --- !
       !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
@@ -1636 +1680 @@
       IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
       ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
-
-      CALL wrk_dealloc( jpi,jpj, zqsr_oce ) 
 #endif
 
@@ -1688 +1730 @@
       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
 
-      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
-      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
+      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
+      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice )
+      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
+      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

@@ -110 +110 @@
       INTEGER  ::   jl                 ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled)
       REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice 
       !!----------------------------------------------------------------------
@@ -197 +196 @@
          ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%]
          !----------------------------------------------------------------------------------------
-         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs )
          CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
 
@@ -203 +202 @@
          CASE( jp_clio )                                       ! CLIO bulk formulation
             ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 
-            ! (zalb_ice) is computed within the bulk routine
-            CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice )
-            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            ! (alb_ice) is computed within the bulk routine
+                                 CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice )
+            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          CASE( jp_core )                                       ! CORE bulk formulation
             ! albedo depends on cloud fraction because of non-linear spectral effects
-            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-            CALL blk_ice_core_flx( t_su, zalb_ice )
-            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+                                 CALL blk_ice_core_flx( t_su, alb_ice )
+            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          CASE ( jp_purecpl )
             ! albedo depends on cloud fraction because of non-linear spectral effects
-            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            ! clem: evap_ice is forced to 0 in coupled mode for now 
-            !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models
-            evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp
-            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          END SELECT
-         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs )
 
          !----------------------------!
@@ -265 +261 @@
       !!----------------------------------------------------------------------
       INTEGER :: ierr
+      INTEGER :: ji, jj
       !!----------------------------------------------------------------------
       IF(lwp) WRITE(numout,*)
@@ -321 +318 @@
       tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
       !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( gphit(ji,jj) > 0._wp ) THEN  ;  rn_amax_2d(ji,jj) = rn_amax_n  ! NH
+            ELSE                             ;  rn_amax_2d(ji,jj) = rn_amax_s  ! SH
+            ENDIF
+        ENDDO
+      ENDDO 
+      !
       nstart = numit  + nn_fsbc      
       nitrun = nitend - nit000 + 1 
@@ -343 +348 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  &
-         &                ln_limdyn, rn_amax, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt  
+         &                ln_limdyn, rn_amax_n, rn_amax_s, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt  
       !!-------------------------------------------------------------------
       !                    
@@ -364 +369 @@
          WRITE(numout,*) '   number of snow layers                                   = ', nlay_s
          WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
-         WRITE(numout,*) '   maximum ice concentration                               = ', rn_amax 
+         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n 
+         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s
          WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
          WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
@@ -579 +585 @@
       sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
       sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
-      sfx_res(:,:) = 0._wp
+      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
       
       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
@@ -595 +601 @@
       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
       hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
-      hfx_err_dif(:,:) = 0._wp   ;
-
+      hfx_err_dif(:,:) = 0._wp
+      wfx_err_sub(:,:) = 0._wp
+      
       afx_tot(:,:) = 0._wp   ;
       afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

@@ -456 +456 @@
       !                                                ! ---------------------------------------- !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-         CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
+         CALL iom_put( "empmr"  , emp    - rnf )                ! upward water flux
+         CALL iom_put( "empbmr" , emp_b  - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
          CALL iom_put( "saltflx", sfx  )                        ! downward salt flux  
                                                                 ! (includes virtual salt flux beneath ice 

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

@@ -68 +68 @@
       !
       rldf = 1     ! For active tracers the 
+      r_fact_lap(:,:,:) = 1.0
 
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
@@ -214 +215 @@
       IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
       IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' )
+      IF( ln_traldf_grif .AND. ln_isfcav         )   &
+           CALL ctl_stop( ' ice shelf and traldf_grif not tested')
       IF( lk_traldf_eiv .AND. .NOT.ln_traldf_iso )   &
            CALL ctl_stop( '          eddy induced velocity on tracers',   &

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

@@ -10 +10 @@
    !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
    !!            3.2  !  2009-04  (G. Madec & NEMO team) 
-   !!            4.0  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
+   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
+   !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
    !!----------------------------------------------------------------------
 
@@ -93 +94 @@
       !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
       !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
+      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
       !!----------------------------------------------------------------------
       !
@@ -101 +103 @@
       REAL(wp) ::   zchl, zcoef, zfact   ! local scalars
       REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
-      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
       REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
+      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
+      REAL(wp) ::   zlogc, zlogc2, zlogc3 
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
-      !!----------------------------------------------------------------------
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d
+      !!--------------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
       !
       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
+      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
       !
       IF( kt == nit000 ) THEN
@@ -183 +186 @@
             !                                             ! ------------------------- !
             ! Set chlorophyl concentration
-            IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume
-               !
-               IF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll
-                  !
-                  CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
-                  !         
-!CDIR COLLAPSE
+            IF( nn_chldta == 1 .OR. nn_chldta == 2 .OR. lk_vvl ) THEN    !*  Variable Chlorophyll or ocean volume
+               !
+               IF( nn_chldta == 1 ) THEN        !*  2D Variable Chlorophyll
+                  !
+                  CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step
+                  DO jk = 1, nksr + 1
+                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 
+                  ENDDO
+                  !
+               ELSE IF( nn_chldta == 2 ) THEN    !*   -3-D Variable Chlorophyll
+                  !
+                  CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step
+!CDIR NOVERRCHK   !
+                  DO jj = 1, jpj
 !CDIR NOVERRCHK
-                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
-!CDIR NOVERRCHK
-                     DO ji = 1, jpi
-                        zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
-                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-                        zekb(ji,jj) = rkrgb(1,irgb)
-                        zekg(ji,jj) = rkrgb(2,irgb)
-                        zekr(ji,jj) = rkrgb(3,irgb)
-                     END DO
-                  END DO
-               ELSE                                            ! Variable ocean volume but constant chrlorophyll
-                  zchl = 0.05                                     ! constant chlorophyll
-                  irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
-                  zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll 
-                  zekg(:,:) = rkrgb(2,irgb)
-                  zekr(:,:) = rkrgb(3,irgb)
+                     DO ji = 1, jpi
+                        zchl    = sf_chl(1)%fnow(ji,jj,1)
+                        zCtot   = 40.6  * zchl**0.459
+                        zze     = 568.2 * zCtot**(-0.746)
+                        IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293)
+                        zlogc   = LOG( zchl )
+                        zlogc2  = zlogc * zlogc
+                        zlogc3  = zlogc * zlogc * zlogc
+                        zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
+                        zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
+                        zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
+                        zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
+                        zCze    = 1.12  * (zchl)**0.803 
+                        DO jk = 1, nksr + 1
+                           zpsi = fsdept(ji,jj,jk) / zze
+                           zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
+                        END DO
+                        !
+                      END DO
+                   END DO
+                     !
+               ELSE                              !* Variable ocean volume but constant chrlorophyll
+                  DO jk = 1, nksr + 1
+                     zchl3d(:,:,jk) = 0.05 
+                  ENDDO
                ENDIF
                !
-               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
+               zcoef  = ( 1. - rn_abs ) / 3.e0                        !  equi-partition in R-G-B
                ze0(:,:,1) = rn_abs  * qsr(:,:)
                ze1(:,:,1) = zcoef * qsr(:,:)
@@ -217 +236 @@
                !
                DO jk = 2, nksr+1
+                  !
+                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of vertical profile of Chl
+!CDIR NOVERRCHK
+                     DO ji = 1, jpi
+                        zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
+                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
+                        zekb(ji,jj) = rkrgb(1,irgb)
+                        zekg(ji,jj) = rkrgb(2,irgb)
+                        zekr(ji,jj) = rkrgb(3,irgb)
+                     END DO
+                  END DO
 !CDIR NOVERRCHK
                   DO jj = 1, jpj
@@ -233 +263 @@
                   END DO
                END DO
-               ! clem: store attenuation coefficient of the first ocean level
-               IF ( ln_qsr_ice ) THEN
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     )
-                        zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) )
-                        zzc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
-                        zzc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
-                        fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2  + zzc3  ) * tmask(ji,jj,2) 
-                     END DO
-                  END DO
-               ENDIF
                !
                DO jk = 1, nksr                                        ! compute and add qsr trend to ta
@@ -251 +269 @@
                zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
                CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution
+               !
+               IF ( ln_qsr_ice ) THEN    ! store attenuation coefficient of the first ocean level
+!CDIR NOVERRCHK
+                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
+!CDIR NOVERRCHK
+                     DO ji = 1, jpi
+                        zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,1) ) )
+                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
+                        zekb(ji,jj) = rkrgb(1,irgb)
+                        zekg(ji,jj) = rkrgb(2,irgb)
+                        zekr(ji,jj) = rkrgb(3,irgb)
+                     END DO
+                  END DO
+                  ! 
+                  DO jj = 1, jpj
+                     DO ji = 1, jpi
+                        zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     )
+                        zc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) )
+                        zc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
+                        zc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
+                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2) 
+                     END DO
+                  END DO
+                  !
+               ENDIF
                !
             ELSE                                                 !*  Constant Chlorophyll
@@ -256 +299 @@
                   qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:)
                END DO
-               ! clem: store attenuation coefficient of the first ocean level
-               IF ( ln_qsr_ice ) THEN
+               !  store attenuation coefficient of the first ocean level
+               IF( ln_qsr_ice ) THEN
                   fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp
                ENDIF
@@ -339 +382 @@
       !
       CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
@@ -405 +448 @@
          WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio
          WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice
-         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta
+         WRITE(numout,*) '      RGB : Chl data (=1/2) or cst value (=0)  nn_chldta  = ', nn_chldta
          WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs
          WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
@@ -429 +472 @@
          IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
          IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2
-         IF( ln_qsr_2bd                      )   nqsr =  3
-         IF( ln_qsr_bio                      )   nqsr =  4
+         IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3
+         IF( ln_qsr_2bd                      )   nqsr =  4
+         IF( ln_qsr_bio                      )   nqsr =  5
          !
          IF(lwp) THEN                   ! Print the choice
             WRITE(numout,*)
             IF( nqsr ==  1 )   WRITE(numout,*) '         R-G-B   light penetration - Constant Chlorophyll'
-            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data '
-            IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration'
-            IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration'
+            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - 2D Chl data '
+            IF( nqsr ==  3 )   WRITE(numout,*) '         R-G-B   light penetration - 3D Chl data '
+            IF( nqsr ==  4 )   WRITE(numout,*) '         2 bands light penetration'
+            IF( nqsr ==  5 )   WRITE(numout,*) '         bio-model light penetration'
          ENDIF
          !
@@ -460 +505 @@
             IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
             !
-            IF( nn_chldta == 1 ) THEN           !* Chl data : set sf_chl structure
+            IF( nn_chldta == 1  .OR. nn_chldta == 2 ) THEN           !* Chl data : set sf_chl structure
                IF(lwp) WRITE(numout,*)
                IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

@@ -177 +177 @@
                   &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
                ! add to the eddy viscosity coef. previously computed
+# if defined key_zdftmx_new
+               ! key_zdftmx_new: New internal wave-driven param: use avs value computed by zdftmx
+               avs (ji,jj,jk) = avs(ji,jj,jk) + zavfs + zavds
+# else
                avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
+# endif
                avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
                avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

@@ -80 +80 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   iikn, iiki, ikt, imkt   ! local integer
-      REAL(wp) ::   zN2_c        ! local scalar
+      INTEGER  ::   ji, jj, jk      ! dummy loop indices
+      INTEGER  ::   iikn, iiki, ikt ! local integer
+      REAL(wp) ::   zN2_c           ! local scalar
       INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
       !!----------------------------------------------------------------------
@@ -117 +117 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               imkt = mikt(ji,jj)
-               IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline 
+               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline 
             END DO
          END DO
@@ -127 +126 @@
             iiki = imld(ji,jj)
             iikn = nmln(ji,jj)
-            imkt = mikt(ji,jj)
-            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Turbocline depth 
-            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Mixed layer depth
-            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
+            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
+            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * ssmask(ji,jj)    ! Mixed layer depth
+            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
          END DO
       END DO
-      IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode
-         CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth
-         CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth
+      ! no need to output in offline mode
+      IF( .NOT.lk_offline ) THEN   
+         IF ( iom_use("mldr10_1") ) THEN
+            IF( ln_isfcav ) THEN
+               CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
+            ELSE
+               CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
+            END IF
+         END IF
+         IF ( iom_use("mldkz5") ) THEN
+            IF( ln_isfcav ) THEN
+               CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
+            ELSE
+               CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
+            END IF
+         END IF
       ENDIF
       

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -357 +357 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcof   = zfact1 * tmask(ji,jj,jk)
+# if defined key_zdftmx_new
+               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
+               zzd_up = zcof * ( MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) )   &  ! upper diagonal
+                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
+               zzd_lw = zcof * ( MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) )   &  ! lower diagonal
+                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+# else
                zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
                   &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
                zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
                   &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+# endif
                   !                                                           ! shear prod. at w-point weightened by mask
                zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
@@ -735 +743 @@
       !
       ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
+# if defined key_zdftmx_new
+      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
+      rn_emin  = 1.e-10_wp
+      rmxl_min = 1.e-03_wp
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+# else
       rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
+# endif
       !
       IF(lwp) THEN                    !* Control print

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

@@ -561 +561 @@
    END SUBROUTINE zdf_tmx_init
 
+#elif defined key_zdftmx_new
+   !!----------------------------------------------------------------------
+   !!   'key_zdftmx_new'               Internal wave-driven vertical mixing
+   !!----------------------------------------------------------------------
+   !!   zdf_tmx       : global     momentum & tracer Kz with wave induced Kz
+   !!   zdf_tmx_init  : global     momentum & tracer Kz with wave induced Kz
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables
+   USE zdf_oce        ! ocean vertical physics variables
+   USE zdfddm         ! ocean vertical physics: double diffusive mixing
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE eosbn2         ! ocean equation of state
+   USE phycst         ! physical constants
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O Manager
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_tmx         ! called in step module 
+   PUBLIC   zdf_tmx_init    ! called in nemogcm module 
+   PUBLIC   zdf_tmx_alloc   ! called in nemogcm module
+
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: wave-driven mixing flag
+
+   !                       !!* Namelist  namzdf_tmx : internal wave-driven mixing *
+   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
+   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency
+   LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F)
+
+   REAL(wp) ::  r1_6 = 1._wp / 6._wp
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_tmx     ! power available from high-mode wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_tmx     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_tmx     ! power available from low-mode, critical slope wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_tmx     ! WKB decay scale for high-mode energy dissipation (m)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_tmx     ! decay scale for low-mode critical slope dissipation (m)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_tmx     ! local energy density available for mixing (W/kg)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_tmx     ! buoyancy flux Kz * N^2 (W/kg)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_tmx    ! vertically integrated buoyancy flux (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity
+
+   !! * Substitutions
+#  include "zdfddm_substitute.h90"
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION zdf_tmx_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_tmx_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE(     ebot_tmx(jpi,jpj),  epyc_tmx(jpi,jpj),  ecri_tmx(jpi,jpj)    ,   &
+      &             hbot_tmx(jpi,jpj),  hcri_tmx(jpi,jpj),  emix_tmx(jpi,jpj,jpk),   &
+      &         bflx_tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk),   & 
+      &         zav_wave(jpi,jpj,jpk), STAT=zdf_tmx_alloc     )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( zdf_tmx_alloc )
+      IF( zdf_tmx_alloc /= 0 )   CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
+   END FUNCTION zdf_tmx_alloc
+
+
+   SUBROUTINE zdf_tmx( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_tmx  ***
+      !!                   
+      !! ** Purpose :   add to the vertical mixing coefficients the effect of
+      !!              breaking internal waves.
+      !!
+      !! ** Method  : - internal wave-driven vertical mixing is given by:
+      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_tmx /( Nu * N^2 )  )
+      !!              where emix_tmx is the 3D space distribution of the wave-breaking 
+      !!              energy and Nu the molecular kinematic viscosity.
+      !!              The function f(Reb) is linear (constant mixing efficiency)
+      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T.
+      !!
+      !!              - Compute emix_tmx, the 3D power density that allows to compute
+      !!              Reb and therefrom the wave-induced vertical diffusivity.
+      !!              This is divided into three components:
+      !!                 1. Bottom-intensified low-mode dissipation at critical slopes
+      !!                     emix_tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx )
+      !!                                   / ( 1. - EXP( - H/hcri_tmx ) ) * hcri_tmx
+      !!              where hcri_tmx is the characteristic length scale of the bottom 
+      !!              intensification, ecri_tmx a map of available power, and H the ocean depth.
+      !!                 2. Pycnocline-intensified low-mode dissipation
+      !!                     emix_tmx(z) = ( epyc_tmx / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
+      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) )
+      !!              where epyc_tmx is a map of available power, and nn_zpyc
+      !!              is the chosen stratification-dependence of the internal wave
+      !!              energy dissipation.
+      !!                 3. WKB-height dependent high mode dissipation
+      !!                     emix_tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx)
+      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_tmx) * e3w(z) )
+      !!              where hbot_tmx is the characteristic length scale of the WKB bottom 
+      !!              intensification, ebot_tmx is a map of available power, and z_wkb is the
+      !!              WKB-stretched height above bottom defined as
+      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) )
+      !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    )
+      !!
+      !!              - update the model vertical eddy viscosity and diffusivity: 
+      !!                     avt  = avt  +    av_wave
+      !!                     avm  = avm  +    av_wave
+      !!                     avmu = avmu + mi(av_wave)
+      !!                     avmv = avmv + mj(av_wave)
+      !!
+      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing:
+      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb)
+      !!
+      !! ** Action  : - Define emix_tmx used to compute internal wave-induced mixing
+      !!              - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing    
+      !!
+      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step 
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   ztpc         ! scalar workspace
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid)
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid)
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('zdf_tmx')
+      !
+      CALL wrk_alloc( jpi,jpj,       zfact, zhdep )
+      CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
+
+      !                          ! ----------------------------- !
+      !                          !  Internal wave-driven mixing  !  (compute zav_wave)
+      !                          ! ----------------------------- !
+      !                             
+      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth,
+      !                                                 using an exponential decay from the seafloor.
+      DO jj = 1, jpj                ! part independent of the level
+         DO ji = 1, jpi
+            zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
+            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_tmx(ji,jj) )  )
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj)
+         END DO
+      END DO
+
+      DO jk = 2, jpkm1              ! complete with the level-dependent part
+         emix_tmx(:,:,jk) = zfact(:,:) * (  EXP( ( fsde3w(:,:,jk  ) - zhdep(:,:) ) / hcri_tmx(:,:) )                      &
+            &                             - EXP( ( fsde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_tmx(:,:) )  ) * wmask(:,:,jk)   &
+            &                          / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) )
+      END DO
+
+      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying 
+      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc
+
+      SELECT CASE ( nn_zpyc )
+
+      CASE ( 1 )               ! Dissipation scales as N (recommended)
+
+         zfact(:,:) = 0._wp
+         DO jk = 2, jpkm1              ! part independent of the level
+            zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         END DO
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
+            END DO
+         END DO
+
+         DO jk = 2, jpkm1              ! complete with the level-dependent part
+            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         END DO
+
+      CASE ( 2 )               ! Dissipation scales as N^2
+
+         zfact(:,:) = 0._wp
+         DO jk = 2, jpkm1              ! part independent of the level
+            zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+         END DO
+
+         DO jj= 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
+            END DO
+         END DO
+
+         DO jk = 2, jpkm1              ! complete with the level-dependent part
+            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+         END DO
+
+      END SELECT
+
+      !                        !* WKB-height dependent mixing: distribute energy over the time-varying 
+      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
+      
+      zwkb(:,:,:) = 0._wp
+      zfact(:,:) = 0._wp
+      DO jk = 2, jpkm1
+         zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         zwkb(:,:,jk) = zfact(:,:)
+      END DO
+
+      DO jk = 2, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
+                                            &           * tmask(ji,jj,jk) / zfact(ji,jj)
+            END DO
+         END DO
+      END DO
+      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1)
+
+      zweight(:,:,:) = 0._wp
+      DO jk = 2, jpkm1
+         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_tmx(:,:) * wmask(:,:,jk)                    &
+            &   * (  EXP( -zwkb(:,:,jk) / hbot_tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) )  )
+      END DO
+
+      zfact(:,:) = 0._wp
+      DO jk = 2, jpkm1              ! part independent of the level
+         zfact(:,:) = zfact(:,:) + zweight(:,:,jk)
+      END DO
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
+         END DO
+      END DO
+
+      DO jk = 2, jpkm1              ! complete with the level-dependent part
+         emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
+            &                                / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) )
+      END DO
+
+
+      ! Calculate molecular kinematic viscosity
+      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  &
+         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0
+      DO jk = 2, jpkm1
+         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
+      END DO
+
+      ! Calculate turbulence intensity parameter Reb
+      DO jk = 2, jpkm1
+         zReb(:,:,jk) = emix_tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
+      END DO
+
+      ! Define internal wave-induced diffusivity
+      DO jk = 2, jpkm1
+         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6
+      END DO
+
+      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
+         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
+                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
+                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s
+         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk)
+      END DO
+
+      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave
+         ztpc = 0._wp
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztpc = ztpc + fse3w(ji,jj,jk) * e1e2t(ji,jj)   &
+                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
+               END DO
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( ztpc )
+         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
+ 
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)'
+            WRITE(numout,*) '~~~~~~~ '
+            WRITE(numout,*)
+            WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW'
+         ENDIF
+      ENDIF
+
+      !                          ! ----------------------- !
+      !                          !   Update  mixing coefs  !                          
+      !                          ! ----------------------- !
+      !      
+      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
+         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  &
+                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   &
+                      &                 ) * wmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "av_ratio", zav_ratio )
+         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing
+            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk)
+            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
+            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
+         END DO
+         !
+      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing
+         DO jk = 2, jpkm1
+            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
+            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
+            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
+         END DO
+      ENDIF
+
+      DO jk = 2, jpkm1              !* update momentum diffusivity at wu and wv points
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1  ! vector opt.
+               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition
+
+      !                             !* output internal wave-driven mixing coefficient
+      CALL iom_put( "av_wave", zav_wave )
+                                    !* output useful diagnostics: N^2, Kz * N^2 (bflx_tmx), 
+                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_tmx), energy density (emix_tmx)
+      IF( iom_use("bflx_tmx") .OR. iom_use("pcmap_tmx") ) THEN
+         bflx_tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)
+         pcmap_tmx(:,:) = 0._wp
+         DO jk = 2, jpkm1
+            pcmap_tmx(:,:) = pcmap_tmx(:,:) + fse3w(:,:,jk) * bflx_tmx(:,:,jk) * wmask(:,:,jk)
+         END DO
+         pcmap_tmx(:,:) = rau0 * pcmap_tmx(:,:)
+         CALL iom_put( "bflx_tmx", bflx_tmx )
+         CALL iom_put( "pcmap_tmx", pcmap_tmx )
+      ENDIF
+      CALL iom_put( "bn2", rn2 )
+      CALL iom_put( "emix_tmx", emix_tmx )
+      
+      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
+
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zdf_tmx')
+      !
+   END SUBROUTINE zdf_tmx
+
+
+   SUBROUTINE zdf_tmx_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_tmx_init  ***
+      !!                     
+      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading
+      !!              of input power maps and decay length scales in netcdf files.
+      !!
+      !! ** Method  : - Read the namzdf_tmx namelist and check the parameters
+      !!
+      !!              - Read the input data in NetCDF files :
+      !!              power available from high-mode wave breaking (mixing_power_bot.nc)
+      !!              power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc)
+      !!              power available from critical slope wave-breaking (mixing_power_cri.nc)
+      !!              WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc)
+      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc)
+      !!
+      !! ** input   : - Namlist namzdf_tmx
+      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc,
+      !!              decay_scale_bot.nc decay_scale_cri.nc
+      !!
+      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
+      !!              - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx
+      !!
+      !! References : de Lavergne et al. 2015, JPO; 2016, in prep.
+      !!         
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   inum         ! local integer
+      INTEGER  ::   ios
+      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars
+      !!
+      NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init')
+      !
+      REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing
+      READ  ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp )
+      !
+      REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing
+      READ  ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namzdf_tmx_new )
+      !
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf_tmx_new : set wave-driven mixing parameters'
+         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc
+         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar
+         WRITE(numout,*) '      Differential internal wave-driven mixing (T) or not (F)   = ', ln_tsdiff
+      ENDIF
+      
+      ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and
+      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should 
+      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6).
+      avmb(:) = 1.4e-6_wp        ! viscous molecular value
+      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_tmx)    
+      avtb_2d(:,:) = 1.e0_wp     ! uniform 
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) '   Force the background value applied to avm & avt in TKE to be everywhere ',   &
+            &               'the viscous molecular value & a very small diffusive value, resp.'
+      ENDIF
+      
+      IF( .NOT.lk_zdfddm )   CALL ctl_stop( 'STOP', 'zdf_tmx_init_new : key_zdftmx_new requires key_zdfddm' )
+      
+      !                             ! allocate tmx arrays
+      IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
+      !
+      !                             ! read necessary fields
+      CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', ebot_tmx, 1 ) 
+      CALL iom_close(inum)
+      !
+      CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', epyc_tmx, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', ecri_tmx, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m]
+      CALL iom_get  (inum, jpdom_data, 'field', hbot_tmx, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m]
+      CALL iom_get  (inum, jpdom_data, 'field', hcri_tmx, 1 )
+      CALL iom_close(inum)
+
+      ebot_tmx(:,:) = ebot_tmx(:,:) * ssmask(:,:)
+      epyc_tmx(:,:) = epyc_tmx(:,:) * ssmask(:,:)
+      ecri_tmx(:,:) = ecri_tmx(:,:) * ssmask(:,:)
+
+      ! Set once for all to zero the first and last vertical levels of appropriate variables
+      emix_tmx (:,:, 1 ) = 0._wp
+      emix_tmx (:,:,jpk) = 0._wp
+      zav_ratio(:,:, 1 ) = 0._wp
+      zav_ratio(:,:,jpk) = 0._wp
+      zav_wave (:,:, 1 ) = 0._wp
+      zav_wave (:,:,jpk) = 0._wp
+
+      zbot = glob_sum( e1e2t(:,:) * ebot_tmx(:,:) )
+      zpyc = glob_sum( e1e2t(:,:) * epyc_tmx(:,:) )
+      zcri = glob_sum( e1e2t(:,:) * ecri_tmx(:,:) )
+      IF(lwp) THEN
+         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Critical slope wave-breaking energy:        ', zcri * 1.e-12_wp, 'TW'
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init')
+      !
+   END SUBROUTINE zdf_tmx_init
+
 #else
    !!----------------------------------------------------------------------

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -338 +338 @@
       !
       IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file
+      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
 
 #if defined key_agrif