Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/DIA

Timestamp:

2012-01-28T17:44:18+01:00 (12 years ago)

Author:

rblod

Message:

Merge of 3.4beta into the trunk

Location:

trunk/NEMOGCM/NEMO/OPA_SRC/DIA

Files:

• diaar5.F90 (modified) (9 diffs)
• diadct.F90 (copied) (copied from branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90) (20 diffs)
• diadimg.F90 (modified) (2 diffs)
• diafwb.F90 (modified) (9 diffs)
• diaharm.F90 (copied) (copied from branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90) (12 diffs)
• diahsb.F90 (modified) (6 diffs)
• diahth.F90 (modified) (11 diffs)
• diaptr.F90 (modified) (21 diffs)
• diawri.F90 (modified) (22 diffs)
• diawri_dimg.h90 (modified) (9 diffs)

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

@@ -19 +19 @@
    USE lib_mpp        ! distribued memory computing library
    USE iom            ! I/O manager library
+   USE timing         ! preformance summary
+   USE wrk_nemo       ! working arrays
 
    IMPLICIT NONE
@@ -65 +67 @@
       !! ** Purpose :   compute and output some AR5 diagnostics
       !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zarea_ssh => wrk_2d_1 , zbotpres => wrk_2d_2   ! 2D workspace
-      USE wrk_nemo, ONLY:   zrhd      => wrk_3d_1 , zrhop    => wrk_3d_2   ! 3D      -
-      USE wrk_nemo, ONLY:   ztsn      => wrk_4d_1                          ! 4D      -
       !
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
@@ -74 +72 @@
       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
       REAL(wp) ::   zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass
+      !
+      REAL(wp), POINTER, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace 
+      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zrhd , zrhop               ! 3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
       !!--------------------------------------------------------------------
-
-      IF( wrk_in_use(2, 1,2) .OR.   &
-          wrk_in_use(3, 1,2) .OR.   &
-          wrk_in_use(4, 1)   ) THEN
-         CALL ctl_stop('dia_ar5: requested workspace arrays unavailable')   ;   RETURN
-      ENDIF
+      IF( nn_timing == 1 )   CALL timing_start('dia_ar5')
+ 
+      CALL wrk_alloc( jpi , jpj              , zarea_ssh , zbotpres )
+      CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    )
+      CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
 
       CALL iom_put( 'cellthc', fse3t(:,:,:) )
@@ -94 +95 @@
       CALL iom_put( 'sshtot', zvolssh / area_tot )
 
-      !                                         ! thermosteric ssh
-      ztsn(:,:,:,jp_tem) = tn (:,:,:)
+      !                     
+      ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
       ztsn(:,:,:,jp_sal) = sn0(:,:,:)
       CALL eos( ztsn, zrhd )                       ! now in situ density using initial salinity
@@ -138 +139 @@
             DO ji = 1, jpi
                zztmp = area(ji,jj) * fse3t(ji,jj,jk)
-               ztemp = ztemp + zztmp * tn(ji,jj,jk)
-               zsal  = zsal  + zztmp * sn(ji,jj,jk)
+               ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
+               zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
             END DO
          END DO
       END DO
       IF( .NOT.lk_vvl ) THEN
-         ztemp = ztemp + SUM( zarea_ssh(:,:) * tn(:,:,1) )
-         zsal  = zsal  + SUM( zarea_ssh(:,:) * sn(:,:,1) )
+         ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
+         zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
       ENDIF
       IF( lk_mpp ) THEN  
@@ -160 +161 @@
       CALL iom_put( 'saltot' , zsal  )
       !
-      IF( wrk_not_released(2, 1,2) .OR.   &
-          wrk_not_released(3, 1,2) .OR.   &
-          wrk_not_released(4, 1)   )   CALL ctl_stop('dia_ar5: failed to release workspace arrays')
+      CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres )
+      CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    )
+      CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5')
       !
    END SUBROUTINE dia_ar5
@@ -173 +176 @@
       !! ** Purpose :   initialization for AR5 diagnostic computation
       !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   wrk_4d_1      ! 4D workspace
-      !
       INTEGER  ::   inum
       INTEGER  ::   ik
@@ -183 +183 @@
       !!----------------------------------------------------------------------
       !
-      IF(wrk_in_use(4, 1) ) THEN
-         CALL ctl_stop('dia_ar5_init: requested workspace array unavailable.')   ;   RETURN
-      ENDIF
-      zsaldta => wrk_4d_1(:,:,:,1:2)
-
+      IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init')
+      !
+      CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
       !                                      ! allocate dia_ar5 arrays
       IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
@@ -221 +219 @@
       ENDIF
       !
-      IF( wrk_not_released(4, 1) )   CALL ctl_stop('dia_ar5_init: failed to release workspace array')
+      CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
       !
    END SUBROUTINE dia_ar5_init

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

@@ -38 +38 @@
   USE dianam          ! build name of file
   USE lib_mpp         ! distributed memory computing library
-#if defined key_ice_lim2 || defined key_ice_lim3
+#if defined key_lim2 || defined key_lim3
   USE ice
 #endif
   USE domvvl
+  USE timing          ! preformance summary
+  USE wrk_nemo        ! working arrays
 
   IMPLICIT NONE
@@ -114 +116 @@
      NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
 
+     IF( nn_timing == 1 )   CALL timing_start('dia_dct_init')
+
      !read namelist
      REWIND( numnam )
@@ -147 +151 @@
      ENDIF
 
-
+     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init')
+     !
   END SUBROUTINE dia_dct_init
  
@@ -161 +166 @@
 
      !! * Local variables
-     INTEGER             :: jsec,            &!loop on sections
-                            iost              !error for opening fileout
-     LOGICAL             :: lldebug =.FALSE.  !debug a section  
-     CHARACTER(len=160)  :: clfileout         !fileout name
+     INTEGER             :: jsec,            &! loop on sections
+                            iost,            &! error for opening fileout
+                            itotal            ! nb_sec_max*nb_type_class*nb_class_max
+     LOGICAL             :: lldebug =.FALSE.  ! debug a section  
+     CHARACTER(len=160)  :: clfileout         ! fileout name
 
      
-     INTEGER , DIMENSION(1):: ish                                       ! tmp array for mpp_sum
-     INTEGER , DIMENSION(3):: ish2                                      !   "
-     REAL(wp), DIMENSION(nb_sec_max*nb_type_class*nb_class_max):: zwork !   "  
-     REAL(wp), DIMENSION(nb_sec_max,nb_type_class,nb_class_max):: zsum  !   "
+     INTEGER , DIMENSION(1)             :: ish   ! tmp array for mpp_sum
+     INTEGER , DIMENSION(3)             :: ish2  !   "
+     REAL(wp), POINTER, DIMENSION(:)    :: zwork !   "  
+     REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   "
 
      !!---------------------------------------------------------------------    
-
+     IF( nn_timing == 1 )   CALL timing_start('dia_dct')
+
+     IF( lk_mpp )THEN
+        itotal = nb_sec_max*nb_type_class*nb_class_max
+        CALL wrk_alloc( itotal                                , zwork ) 
+        CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum  )
+     ENDIF    
+ 
      IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
          WRITE(numout,*) " "
@@ -189 +202 @@
            !debug this section computing ?
            lldebug=.FALSE.
-!           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. 
-           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 ) lldebug=.TRUE. 
+           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. 
 
            !Compute transport through section  
@@ -226 +238 @@
      ENDIF
 
+     IF( lk_mpp )THEN
+        itotal = nb_sec_max*nb_type_class*nb_class_max
+        CALL wrk_dealloc( itotal                                , zwork ) 
+        CALL wrk_dealloc( nb_sec_max,nb_type_class,nb_class_max , zsum  )
+     ENDIF    
+
+     IF( nn_timing == 1 )   CALL timing_stop('dia_dct')
+     !
   END SUBROUTINE dia_dct
 
@@ -250 +270 @@
      TYPE(POINT_SECTION),DIMENSION(nb_point_max)  ::coordtemp !contains listpoints coordinates 
                                                               !read in the file
-     INTEGER,DIMENSION(nb_point_max)  ::directemp             !contains listpoints directions
+     INTEGER, POINTER, DIMENSION(:) :: directemp              !contains listpoints directions
                                                               !read in the files
      LOGICAL :: llbon                                       ,&!local logical
                 lldebug                                       !debug the section
      !!-------------------------------------------------------------------------------------
+     CALL wrk_alloc( nb_point_max, directemp )
 
      !open input file
@@ -381 +402 @@
            ENDIF
 
+              IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
+              WRITE(narea+200,*)'avant secs(jsec)%nb_point iptloc ',secs(jsec)%nb_point,iptloc
+              DO jpt = 1,iptloc
+                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
+                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
+                 WRITE(narea+200,*)'avant # I J : ',iiglo,ijglo
+              ENDDO
+              ENDIF
+
            !remove redundant points between processors
            !------------------------------------------
@@ -390 +420 @@
               CALL removepoints(secs(jsec),'J','bot_list',lldebug)
            ENDIF
+           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
+              WRITE(narea+200,*)'apres secs(jsec)%nb_point iptloc ',secs(jsec)%nb_point,iptloc
+              DO jpt = 1,secs(jsec)%nb_point
+                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
+                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
+                 WRITE(narea+200,*)'apres # I J : ',iiglo,ijglo
+              ENDDO
+           ENDIF
 
            !debug
@@ -395 +433 @@
            IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
               WRITE(numout,*)"      List of points after removepoints:"
+              iptloc = secs(jsec)%nb_point
               DO jpt = 1,iptloc
                  iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
@@ -411 +450 @@
      nb_sec = jsec-1   !number of section read in the file
 
+     CALL wrk_dealloc( nb_point_max, directemp )
+     !
   END SUBROUTINE readsec
 
@@ -436 +477 @@
                                  ! isgn=-1 : scan listpoint from end to start 
                 istart,iend      !first and last points selected in listpoint
-     INTEGER :: jpoint   =0      !loop on list points
-     INTEGER,DIMENSION(nb_point_max)   :: idirec !contains temporary sec%direction
-     INTEGER,DIMENSION(2,nb_point_max) :: icoord !contains temporary sec%listpoint
+     INTEGER :: jpoint           !loop on list points
+     INTEGER, POINTER, DIMENSION(:)   :: idirec !contains temporary sec%direction
+     INTEGER, POINTER, DIMENSION(:,:) :: icoord !contains temporary sec%listpoint
      !----------------------------------------------------------------------------
+     CALL wrk_alloc(    nb_point_max, idirec )
+     CALL wrk_alloc( 2, nb_point_max, icoord )
+
      IF( ld_debug )WRITE(numout,*)'      -------------------------'
      IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
@@ -467 +511 @@
      sec%direction            = 0
 
-
      jpoint=iextr+isgn
-     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point  .AND. &
-        icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )
-        jpoint=jpoint+isgn
-     ENDDO
+     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
+         IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
+         ELSE                                                                               ; EXIT
+         ENDIF
+     ENDDO 
 
      IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
      ELSE                        ; istart=1        ; iend=jpoint+1
      ENDIF
+
      sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
      sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
@@ -487 +532 @@
      ENDIF
 
+     CALL wrk_dealloc(    nb_point_max, idirec )
+     CALL wrk_dealloc( 2, nb_point_max, icoord )
   END SUBROUTINE removepoints
 
@@ -536 +583 @@
 
      TYPE(POINT_SECTION) :: k
-     REAL(wp),DIMENSION(nb_type_class,nb_class_max)::zsum
+     REAL(wp), POINTER, DIMENSION(:,:):: zsum ! 2D work array
      !!--------------------------------------------------------
+     CALL wrk_alloc( nb_type_class , nb_class_max , zsum   )
 
      IF( ld_debug )WRITE(numout,*)'      Compute transport'
@@ -746 +794 @@
            ENDDO !end of loop on the density classes
 
-#if defined key_ice_lim
+#if defined key_lim2 || defined key_lim3
 
            !ICE CASE    
@@ -812 +860 @@
         sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6
         IF( sec%llstrpond ) THEN
-           IF( zsum(1,jclass) .NE. 0 ) THEN
-              sec%transport(3,jclass)=sec%transport(3,jclass)+zsum(3,jclass)/zsum(1,jclass)
-              sec%transport(5,jclass)=sec%transport(5,jclass)+zsum(5,jclass)/zsum(1,jclass)
-              sec%transport(7,jclass)=sec%transport(7,jclass)+zsum(7,jclass)
-              sec%transport(9,jclass)=sec%transport(9,jclass)+zsum(9,jclass)
-           ELSE
-              sec%transport(3,jclass)=0.
-              sec%transport(5,jclass)=0.
-              sec%transport(7,jclass)=0.
-              sec%transport(9,jclass)=0.
+           IF( zsum(1,jclass) .NE. 0._wp ) THEN
+              sec%transport( 3,jclass) = sec%transport( 3,jclass) + zsum( 3,jclass)/zsum(1,jclass)
+              sec%transport( 5,jclass) = sec%transport( 5,jclass) + zsum( 5,jclass)/zsum(1,jclass)
+              sec%transport( 7,jclass) = sec%transport( 7,jclass) + zsum( 7,jclass)
+              sec%transport( 9,jclass) = sec%transport( 9,jclass) + zsum( 9,jclass)
            ENDIF
-           IF( zsum(2,jclass) .NE. 0 )THEN
-              sec%transport( 4,jclass)=sec%transport( 4,jclass)+zsum( 4,jclass)/zsum(2,jclass)
-              sec%transport( 6,jclass)=sec%transport( 6,jclass)+zsum( 6,jclass)/zsum(2,jclass)
-              sec%transport( 8,jclass)=sec%transport( 8,jclass)+zsum( 8,jclass)
-              sec%transport(10,jclass)=sec%transport(10,jclass)+zsum(10,jclass)
-           ELSE
-              sec%transport( 4,jclass)=0.
-              sec%transport( 6,jclass)=0.
-              sec%transport( 8,jclass)=0.
-              sec%transport(10,jclass)=0.
+           IF( zsum(2,jclass) .NE. 0._wp )THEN
+              sec%transport( 4,jclass) = sec%transport( 4,jclass) + zsum( 4,jclass)/zsum(2,jclass)
+              sec%transport( 6,jclass) = sec%transport( 6,jclass) + zsum( 6,jclass)/zsum(2,jclass)
+              sec%transport( 8,jclass) = sec%transport( 8,jclass) + zsum( 8,jclass)
+              sec%transport(10,jclass) = sec%transport(10,jclass) + zsum(10,jclass)
            ENDIF
         ELSE
-           sec%transport( 3,jclass)=0.
-           sec%transport( 4,jclass)=0.
-           sec%transport( 5,jclass)=0.
-           sec%transport( 6,jclass)=0.
-           sec%transport( 7,jclass)=0.
-           sec%transport( 8,jclass)=0.
-           sec%transport(10,jclass)=0.
+           sec%transport( 3,jclass) = 0._wp
+           sec%transport( 4,jclass) = 0._wp
+           sec%transport( 5,jclass) = 0._wp
+           sec%transport( 6,jclass) = 0._wp
+           sec%transport( 7,jclass) = 0._wp
+           sec%transport( 8,jclass) = 0._wp
+           sec%transport(10,jclass) = 0._wp
         ENDIF
      ENDDO   
@@ -852 +890 @@
      ENDIF
 
+     CALL wrk_dealloc( nb_type_class , nb_class_max , zsum   )
+     !
   END SUBROUTINE transport
   
@@ -872 +912 @@
      !!------------------------------------------------------------- 
      !!arguments
-     INTEGER, INTENT(IN)          :: kt         ! time-step
-     TYPE(SECTION), INTENT(INOUT) :: sec        ! section to write   
-     INTEGER ,INTENT(IN)          :: ksec       ! section number
+     INTEGER, INTENT(IN)          :: kt          ! time-step
+     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
+     INTEGER ,INTENT(IN)          :: ksec        ! section number
 
      !!local declarations
-     REAL(wp) ,DIMENSION(nb_type_class):: zsumclass
-     INTEGER               :: jcl,ji            ! Dummy loop
-     CHARACTER(len=2)      :: classe            ! Classname 
-     REAL(wp)              :: zbnd1,zbnd2       ! Class bounds
-     REAL(wp)              :: zslope            ! section's slope coeff
+     INTEGER               :: jcl,ji             ! Dummy loop
+     CHARACTER(len=2)      :: classe             ! Classname 
+     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
+     REAL(wp)              :: zslope             ! section's slope coeff
+     !
+     REAL(wp), POINTER, DIMENSION(:):: zsumclass ! 1D workspace 
      !!------------------------------------------------------------- 
-      
+     CALL wrk_alloc(nb_type_class , zsumclass )  
+
      zsumclass(:)=0._wp
      zslope = sec%slopeSection       
@@ -996 +1038 @@
 119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
 
+     CALL wrk_dealloc(nb_type_class , zsumclass )  
   END SUBROUTINE dia_dct_wri
 

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90

@@ -10 +10 @@
    USE in_out_manager  ! I/O manager
    USE daymod          ! calendar
+   USE lib_mpp
 
    IMPLICIT NONE
@@ -40 +41 @@
       !!---------------------------------------------------------------------
       !
-      ALLOCATE( z42d(jpi,jpj), z4dep(jpk), STAT=dia_wri_dimg_alloc )
-      !
-      IF( lk_mpp                  )   CALL mpp_sum ( dia_wri_dimg_alloc )
-      IF( dia_wri_dimg_alloc /= 0 )   CALL ctl_warn('dia_wri_dimg_alloc: allocation of array failed.')
+      IF( .NOT. ALLOCATED( z42d ) )THEN
+
+         ALLOCATE( z42d(jpi,jpj), z4dep(jpk), STAT=dia_wri_dimg_alloc )
+
+         IF( lk_mpp                  )   CALL mpp_sum ( dia_wri_dimg_alloc )
+         IF( dia_wri_dimg_alloc /= 0 )   CALL ctl_warn('dia_wri_dimg_alloc: allocation of array failed.')
+
+      ELSE
+
+         dia_wri_dimg_alloc = 0
+
+      ENDIF
       !
   END FUNCTION dia_wri_dimg_alloc

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90

@@ -22 +22 @@
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distributed memory computing library
+   USE timing          ! preformance summary
 
    IMPLICIT NONE
@@ -61 +62 @@
       REAL(wp) ::  zsm0, zfwfnew
       !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )   CALL timing_start('dia_fwb')
 
       ! Mean global salinity
@@ -80 +82 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-                  a_salb = a_salb + ( sb(ji,jj,jk) - zsm0 ) * zwei
+                  a_salb = a_salb + ( tsb(ji,jj,jk,jp_sal) - zsm0 ) * zwei
                END DO
             END DO
@@ -106 +108 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zwei  = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-                  a_saln = a_saln + ( sn(ji,jj,jk) - zsm0 ) * zwei
+                  a_saln = a_saln + ( tsn(ji,jj,jk,jp_sal) - zsm0 ) * zwei
                   zvol  = zvol  + zwei
                END DO
@@ -174 +176 @@
          END SELECT
          ! 
-         DO ji = mi0(ii0), mi1(ii1)
+         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
             DO jj = mj0(ij0), mj1(ij1)
                DO jk = 1, jpk 
-                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
-                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
-                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
+                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
 
                   IF( un(ji,jj,jk) > 0.e0 ) THEN 
@@ -221 +223 @@
          END SELECT
          ! 
-         DO ji = mi0(ii0), mi1(ii1)
+         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
             DO jj = mj0(ij0), mj1(ij1)
                DO jk = 1, jpk 
-                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
-                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
-                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
+                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
                   
                   IF( un(ji,jj,jk) > 0.e0 ) THEN 
@@ -268 +270 @@
          END SELECT
          ! 
-         DO ji = mi0(ii0), mi1(ii1)
+         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
             DO jj = mj0(ij0), mj1(ij1)
                DO jk = 1, jpk 
-                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
-                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
-                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
+                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
                   
                   IF( un(ji,jj,jk) > 0.e0 ) THEN 
@@ -315 +317 @@
          END SELECT
          ! 
-         DO ji = mi0(ii0), mi1(ii1)
+         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
             DO jj = mj0(ij0), mj1(ij1)
                DO jk = 1, jpk
-                  zt = 0.5 * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
-                  zs = 0.5 * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) )
-                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
+                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+                  zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
                   
                   IF( un(ji,jj,jk) > 0.e0 ) THEN 
@@ -438 +440 @@
       ENDIF
 
+      IF( nn_timing == 1 )   CALL timing_start('dia_fwb')
+
  9005 FORMAT(1X,A,ES24.16)
  9010 FORMAT(1X,A,ES12.5,A,F10.5,A)

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

@@ -16 +16 @@
    USE dynspg_oce
    USE dynspg_ts
-   USE surdetermine
    USE daymod
    USE tide_mod
    USE iom 
+   USE timing          ! preformance summary
+   USE wrk_nemo        ! working arrays
 
    IMPLICIT NONE
    PRIVATE
+
+   LOGICAL, PUBLIC, PARAMETER :: lk_diaharm  = .TRUE.
    
-   INTEGER, PARAMETER :: nb_harmo_max=9
-
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm  = .TRUE.
-
-   INTEGER ::      & !! namelist variables
-       nit000_han=1, & ! First time step used for harmonic analysis
-       nitend_han=1, & ! Last time step used for harmonic analysis
-       nstep_han=1,  & ! Time step frequency for harmonic analysis
-       nb_ana          ! Number of harmonics to analyse
-
-   CHARACTER (LEN=4), DIMENSION(nb_harmo_max) ::   &
+   INTEGER, PARAMETER :: jpincomax    = 2.*jpmax_harmo
+   INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24
+
+   INTEGER ::                            & !! namelist variables
+                         nit000_han = 1, & ! First time step used for harmonic analysis
+                         nitend_han = 1, & ! Last time step used for harmonic analysis
+                         nstep_han  = 1, & ! Time step frequency for harmonic analysis
+                         nb_ana            ! Number of harmonics to analyse
+
+   INTEGER , ALLOCATABLE, DIMENSION(:)       :: name
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
+   REAL(wp), ALLOCATABLE, DIMENSION(:)       :: ana_freq, vt, ut, ft
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: out_eta, &
+                                                out_u  , &
+                                                out_v
+
+   INTEGER :: ninco, nsparse
+   INTEGER ,       DIMENSION(jpdimsparse)         :: njsparse, nisparse
+   INTEGER , SAVE, DIMENSION(jpincomax)           :: ipos1
+   REAL(wp),       DIMENSION(jpdimsparse)         :: valuesparse
+   REAL(wp),       DIMENSION(jpincomax)           :: ztmp4 , ztmp7
+   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
+   REAL(wp), SAVE, DIMENSION(jpincomax)           :: zpivot
+
+   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   &
        tname         ! Names of tidal constituents ('M2', 'K1',...)
 
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
-   REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, vt, ut, ft
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta, &
-              out_u, &
-                 out_v
-   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: name
 
 !! * Routine accessibility
@@ -67 +78 @@
       !! * Local declarations 
       INTEGER :: jh, nhan, jk, ji
-      NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, nb_ana, tname
-
-      !!----------------------------------------------------------------------
-
-      ! Read namelist parameters:
-      ! -------------------------
+
+      NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
+         WRITE(numout,*) '~~~~~~~ '
+      ENDIF
+      !
+      CALL tide_init_Wave
+      !
+      tname(:)=''
+      !
+      ! Read Namelist nam_diaharm
       REWIND ( numnam )
       READ   ( numnam, nam_diaharm )
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
-
-      IF(lwp) WRITE(numout,*) 'First time step used for analysis:  nit000_han= ', nit000_han
-      IF(lwp) WRITE(numout,*) 'Last  time step used for analysis:  nitend_han= ', nitend_han
-      IF(lwp) WRITE(numout,*) 'Time step frequency for harmonic analysis:  nstep_han= ', nstep_han
-
-      IF (nb_ana > nb_harmo_max) THEN
-        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
-                                & nb_ana must be lower than nb_harmo_max, stop'
-        IF(lwp) WRITE(numout,*) 'nb_harmo_max= ', nb_harmo_max
-        nstop = nstop + 1
+      !
+      IF(lwp) THEN
+         WRITE(numout,*) 'First time step used for analysis:  nit000_han= ', nit000_han
+         WRITE(numout,*) 'Last  time step used for analysis:  nitend_han= ', nitend_han
+         WRITE(numout,*) 'Time step frequency for harmonic analysis:  nstep_han= ', nstep_han
       ENDIF
 
@@ -94 +105 @@
       ! ----------------------------------------------
       IF (nit000 > nit000_han) THEN
-        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
-                                & nit000_han must be greater than nit000, stop'
-        IF(lwp) WRITE(numout,*) 'restart capability not implemented'
+        IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nit000_han must be greater than nit000, stop'
+        IF(lwp) WRITE(numout,*) ' restart capability not implemented'
         nstop = nstop + 1
       ENDIF
       IF (nitend < nitend_han) THEN
-        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
-                                & nitend_han must be lower than nitend, stop'
-        IF(lwp) WRITE(numout,*) 'restart capability not implemented'
+        IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nitend_han must be lower than nitend, stop'
+        IF(lwp) WRITE(numout,*) ' restart capability not implemented'
         nstop = nstop + 1
       ENDIF
 
       IF (MOD(nitend_han-nit000_han+1,nstep_han).NE.0) THEN
-        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
-                                & analysis time span must be a multiple of nstep_han, stop'
+        IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : analysis time span must be a multiple of nstep_han, stop'
         nstop = nstop + 1
       END IF
 
-      CALL tide_init_Wave
+      nb_ana=0
+      DO jk=1,jpmax_harmo
+         DO ji=1,jpmax_harmo
+            IF(TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN
+               nb_ana=nb_ana+1
+            ENDIF
+         END DO
+      ENDDO
+      !
+      IF(lwp) THEN
+         WRITE(numout,*) '        Namelist nam_diaharm'
+         WRITE(numout,*) '        nb_ana    = ', nb_ana
+         CALL flush(numout)
+      ENDIF
+      !
+      IF (nb_ana > jpmax_harmo) THEN
+        IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nb_ana must be lower than jpmax_harmo, stop'
+        IF(lwp) WRITE(numout,*) ' jpmax_harmo= ', jpmax_harmo
+        nstop = nstop + 1
+      ENDIF
 
       ALLOCATE(name    (nb_ana))
@@ -162 +189 @@
 
       !! * Local declarations
-      INTEGER :: ji, jj, jh, jc, nhc
+      INTEGER  :: ji, jj, jh, jc, nhc
       REAL(wp) :: ztime, ztemp
+      !!--------------------------------------------------------------------
+      IF( nn_timing == 1 )   CALL timing_start('dia_harm')
 
       IF ( kt .EQ. nit000 ) CALL dia_harm_init
@@ -170 +199 @@
            (MOD(kt,nstep_han).EQ.0) ) THEN
 
-        ztime = kt*rdt 
+        ztime = (kt-nit000+1)*rdt 
        
         nhc = 0
@@ -202 +231 @@
       IF ( kt .EQ. nitend_han ) CALL dia_harm_end
 
+      IF( nn_timing == 1 )   CALL timing_stop('dia_harm')
  
    END SUBROUTINE dia_harm
@@ -223 +253 @@
       REAL(wp) :: ztime, ztime_ini, ztime_end
       REAL(wp) :: X1,X2
-      REAL(wp), DIMENSION(jpi,jpj,nb_harmo_max,2) :: ana_amp
-
-
-        IF(lwp) WRITE(numout,*)
-        IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
-        IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-
-        ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis
-        ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis
-        nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
-
-        ninco = 2*nb_ana
-
-        ksp = 0
-        keq = 0        
-        DO jn = 1, nhan
-          ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1)
-          keq = keq + 1
-          kun = 0
-          DO jh = 1,nb_ana
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ana_amp
+      !!--------------------------------------------------------------------
+      CALL wrk_alloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+
+      ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis
+      ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis
+      nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
+
+      ninco = 2*nb_ana
+
+      ksp = 0
+      keq = 0        
+      DO jn = 1, nhan
+         ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1)
+         keq = keq + 1
+         kun = 0
+         DO jh = 1,nb_ana
             DO jc = 1,2
-              kun = kun + 1
-              ksp = ksp + 1
-              nisparse(ksp) = keq
-              njsparse(ksp) = kun
-              valuesparse(ksp)= &
-                 +(     MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
+               kun = kun + 1
+               ksp = ksp + 1
+               nisparse(ksp) = keq
+               njsparse(ksp) = kun
+               valuesparse(ksp)= &
+                   +(   MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
                    +(1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
             END DO
-          END DO
-        END DO
-
-        nsparse=ksp
-
-        ! Elevation:
-        DO jj = 1, jpj
-          DO ji = 1, jpi
+         END DO
+      END DO
+
+      nsparse=ksp
+
+      ! Elevation:
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             ! Fill input array
             kun=0
             DO jh = 1,nb_ana
-              DO jc = 1,2
-                kun = kun + 1
-                tmp4(kun)=ana_temp(ji,jj,kun,1)
-              ENDDO
+               DO jc = 1,2
+                  kun = kun + 1
+                  ztmp4(kun)=ana_temp(ji,jj,kun,1)
+               ENDDO
             ENDDO
 
@@ -273 +304 @@
             ! Fill output array
             DO jh = 1, nb_ana
-              ana_amp(ji,jj,jh,1)=tmp7((jh-1)*2+1)
-              ana_amp(ji,jj,jh,2)=tmp7((jh-1)*2+2)
+               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
+               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
             END DO
-          END DO
-        END DO
-
-        ALLOCATE(out_eta(jpi,jpj,2*nb_ana))
-        ALLOCATE(out_u  (jpi,jpj,2*nb_ana))
-        ALLOCATE(out_v  (jpi,jpj,2*nb_ana))
-
-
-        DO jj = 1, jpj
-          DO ji = 1, jpi
+         END DO
+      END DO
+
+      ALLOCATE(out_eta(jpi,jpj,2*nb_ana))
+      ALLOCATE(out_u  (jpi,jpj,2*nb_ana))
+      ALLOCATE(out_v  (jpi,jpj,2*nb_ana))
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             DO jh = 1, nb_ana 
-                X1=ana_amp(ji,jj,jh,1)
-                X2=-ana_amp(ji,jj,jh,2)
-                out_eta(ji,jj,jh)=X1 * tmask(ji,jj,1)
-                out_eta(ji,jj,nb_ana+jh)=X2 * tmask(ji,jj,1)
+               X1=ana_amp(ji,jj,jh,1)
+               X2=-ana_amp(ji,jj,jh,2)
+               out_eta(ji,jj,jh)=X1 * tmask(ji,jj,1)
+               out_eta(ji,jj,nb_ana+jh)=X2 * tmask(ji,jj,1)
             ENDDO
-          ENDDO
-        ENDDO
-
-        ! ubar:
-        DO jj = 1, jpj
-          DO ji = 1, jpi
+         ENDDO
+      ENDDO
+
+      ! ubar:
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             ! Fill input array
             kun=0
             DO jh = 1,nb_ana
-              DO jc = 1,2
-                kun = kun + 1
-                tmp4(kun)=ana_temp(ji,jj,kun,2)
-              ENDDO
+               DO jc = 1,2
+                  kun = kun + 1
+                  ztmp4(kun)=ana_temp(ji,jj,kun,2)
+               ENDDO
             ENDDO
 
@@ -311 +341 @@
             ! Fill output array
             DO jh = 1, nb_ana
-              ana_amp(ji,jj,jh,1)=tmp7((jh-1)*2+1)
-              ana_amp(ji,jj,jh,2)=tmp7((jh-1)*2+2)
+               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
+               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
             END DO
 
-          END DO
-        END DO
-
-        DO jj = 1, jpj
-          DO ji = 1, jpi
+         END DO
+      END DO
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             DO jh = 1, nb_ana 
-                X1=ana_amp(ji,jj,jh,1)
-                X2=-ana_amp(ji,jj,jh,2)
-                 out_u(ji,jj,jh) = X1 * umask(ji,jj,1)
-                 out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1)
+               X1=ana_amp(ji,jj,jh,1)
+               X2=-ana_amp(ji,jj,jh,2)
+               out_u(ji,jj,jh) = X1 * umask(ji,jj,1)
+               out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1)
             ENDDO
-          ENDDO
-        ENDDO
-
-        ! vbar:
-        DO jj = 1, jpj
-          DO ji = 1, jpi
-              ! Fill input array
-              kun=0
-              DO jh = 1,nb_ana
-                DO jc = 1,2
+         ENDDO
+      ENDDO
+
+      ! vbar:
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ! Fill input array
+            kun=0
+            DO jh = 1,nb_ana
+               DO jc = 1,2
                   kun = kun + 1
-                  tmp4(kun)=ana_temp(ji,jj,kun,3)
-                ENDDO
-              ENDDO
-
-              CALL SUR_DETERMINE(jj+1)
-
-              ! Fill output array
-              DO jh = 1, nb_ana
-                ana_amp(ji,jj,jh,1)=tmp7((jh-1)*2+1)
-                ana_amp(ji,jj,jh,2)=tmp7((jh-1)*2+2)
-              END DO
-
-          END DO
-        END DO
-
-        DO jj = 1, jpj
-          DO ji = 1, jpi
+                  ztmp4(kun)=ana_temp(ji,jj,kun,3)
+               ENDDO
+            ENDDO
+
+            CALL SUR_DETERMINE(jj+1)
+
+            ! Fill output array
+            DO jh = 1, nb_ana
+               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
+               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
+            END DO
+
+         END DO
+      END DO
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             DO jh = 1, nb_ana 
-                X1=ana_amp(ji,jj,jh,1)
-                X2=-ana_amp(ji,jj,jh,2)
-                 out_v(ji,jj,jh)=X1 * vmask(ji,jj,1)
-                 out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1)
+               X1=ana_amp(ji,jj,jh,1)
+               X2=-ana_amp(ji,jj,jh,2)
+               out_v(ji,jj,jh)=X1 * vmask(ji,jj,1)
+               out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1)
             ENDDO
-          ENDDO
-        ENDDO
-
-        CALL dia_wri_harm ! Write results in files
-
-  END SUBROUTINE dia_harm_end
-
-  SUBROUTINE dia_wri_harm
+         ENDDO
+      ENDDO
+
+      CALL dia_wri_harm ! Write results in files
+      CALL wrk_dealloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
+      !
+   END SUBROUTINE dia_harm_end
+
+   SUBROUTINE dia_wri_harm
       !!--------------------------------------------------------------------
       !!                 ***  ROUTINE dia_wri_harm  ***
@@ -382 +413 @@
       CHARACTER(LEN=lc) :: cltext
       CHARACTER(LEN=lc) ::   &
-         cdfile_name_T    ,   & ! name of the file created (T-points)
-         cdfile_name_U    ,   & ! name of the file created (U-points)
-         cdfile_name_V          ! name of the file created (V-points)
-      INTEGER  ::   jh       
+         cdfile_name_T   ,   & ! name of the file created (T-points)
+         cdfile_name_U   ,   & ! name of the file created (U-points)
+         cdfile_name_V         ! name of the file created (V-points)
+      INTEGER  ::   jh
       !!----------------------------------------------------------------------
 
@@ -395 +426 @@
 
       IF(lwp) WRITE(numout,*) '  '
-      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 
+      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
 #if defined key_dimgout
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T)
@@ -442 +473 @@
 #endif
 
-  END SUBROUTINE dia_wri_harm
+   END SUBROUTINE dia_wri_harm
+
+   SUBROUTINE SUR_DETERMINE(init)
+   !!---------------------------------------------------------------------------------
+   !!                      *** ROUTINE SUR_DETERMINE ***
+   !!    
+   !!    
+   !!       
+   !!---------------------------------------------------------------------------------
+   INTEGER, INTENT(in) :: init 
+               
+   INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
+   REAL(wp)                        :: zval1, zval2, zx1
+   REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
+   INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
+   !---------------------------------------------------------------------------------
+   CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
+   CALL wrk_alloc( jpincomax , ipos2 , ipivot        )
+            
+   IF( init==1 )THEN
+
+      IF( nsparse .GT. jpdimsparse ) &
+         CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
+
+      IF( ninco .GT. jpincomax ) &
+         CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
+
+      ztmp3(:,:)=0.e0
+
+      DO jk1_sd = 1, nsparse
+         DO jk2_sd = 1, nsparse
+
+            nisparse(jk2_sd)=nisparse(jk2_sd)
+            njsparse(jk2_sd)=njsparse(jk2_sd)
+
+            IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
+               ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
+                                                        + valuesparse(jk1_sd)*valuesparse(jk2_sd)
+            ENDIF
+
+         ENDDO
+      ENDDO
+
+      DO jj_sd = 1 ,ninco
+          ipos1(jj_sd) = jj_sd
+          ipos2(jj_sd) = jj_sd
+      ENDDO
+
+      DO ji_sd = 1 , ninco
+
+         !find greatest non-zero pivot:
+         zval1 = ABS(ztmp3(ji_sd,ji_sd))
+
+         ipivot(ji_sd) = ji_sd
+         DO jj_sd = ji_sd, ninco
+            zval2 = ABS(ztmp3(ji_sd,jj_sd))
+            IF( zval2.GE.zval1 )THEN
+               ipivot(ji_sd) = jj_sd
+               zval1         = zval2
+            ENDIF
+         ENDDO
+
+         DO ji1_sd = 1, ninco
+            zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd)
+            zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd))
+            ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd)
+            ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
+         ENDDO
+
+         ipos2(ji_sd)         = ipos1(ipivot(ji_sd))
+         ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
+         ipos1(ji_sd)         = ipos2(ji_sd)
+         ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
+         zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd)
+         DO jj_sd = 1, ninco
+            ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
+         ENDDO
+
+         DO ji2_sd = ji_sd+1, ninco
+            zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
+            DO jj_sd=1,ninco
+               ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
+            ENDDO
+         ENDDO
+
+      ENDDO
+
+   ENDIF ! End init==1
+
+   DO ji_sd = 1, ninco
+      ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
+      DO ji2_sd = ji_sd+1, ninco
+         ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
+      ENDDO
+   ENDDO
+
+   !system solving: 
+   ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
+   ji_sd = ninco
+   DO ji_sd = ninco-1, 1, -1
+      zx1=0.
+      DO jj_sd = ji_sd+1, ninco
+         zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
+      ENDDO
+      ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
+   ENDDO
+
+   DO jj_sd =1, ninco
+      ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
+   ENDDO
+
+
+   CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
+   CALL wrk_dealloc( jpincomax , ipos2 , ipivot        )
+
+  END SUBROUTINE SUR_DETERMINE
+
 
 #else

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

@@ -18 +18 @@
    USE lib_mpp         ! distributed memory computing library
    USE trabbc          ! bottom boundary condition
+   USE obc_par         ! (for lk_obc)
    USE bdy_par         ! (for lk_bdy)
-   USE obc_par         ! (for lk_obc)
+   USE timing          ! preformance summary
 
    IMPLICIT NONE
@@ -72 +73 @@
       REAL(dp)   ::   z_frc_trd_v                 !    -     -
       !!---------------------------------------------------------------------------
+      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')
 
       ! ------------------------- !
@@ -107 +109 @@
          ! heat content variation
          zdiff_hc = zdiff_hc + SUM( surf(:,:) * tmask(:,:,jk)          &
-            &                       * ( fse3t_n(:,:,jk) * tn(:,:,jk)   &
+            &                       * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   &
             &                           - hc_loc_ini(:,:,jk) ) )
          ! salt content variation
          zdiff_sc = zdiff_sc + SUM( surf(:,:) * tmask(:,:,jk)          &
-            &                       * ( fse3t_n(:,:,jk) * sn(:,:,jk)   &
+            &                       * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   &
             &                           - sc_loc_ini(:,:,jk) ) )
       ENDDO
@@ -139 +141 @@
       IF ( kt == nitend ) CLOSE( numhsb )
 
+      IF( nn_timing == 1 )   CALL timing_stop('dia_hsb')
+
 9020  FORMAT(I5,11D15.7)
       !
@@ -205 +209 @@
       WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
       WRITE(numout,*) "~~~~~~~  output written in the 'heat_salt_volume_budgets.txt' ASCII file"
-      IF( lk_obc .OR. lk_bdy) THEN
+      IF( lk_obc .or. lk_bdy ) THEN
          CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
       ENDIF
@@ -248 +252 @@
       ! 4 - initial conservation variables !
       ! ---------------------------------- !
-      ssh_ini(:,:) = sshn(:,:)                               ! initial ssh
+      ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh
       DO jk = 1, jpk
-         e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                ! initial vertical scale factors
-         hc_loc_ini(:,:,jk) = tn(:,:,jk) * fse3t_n(:,:,jk)   ! initial heat content
-         sc_loc_ini(:,:,jk) = sn(:,:,jk) * fse3t_n(:,:,jk)   ! initial salt content
+         e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors
+         hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content
+         sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
       END DO
       frc_v = 0.d0                                           ! volume       trend due to forcing

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90

@@ -23 +23 @@
    USE lib_mpp         ! MPP library
    USE iom             ! I/O library
+   USE timing          ! preformance summary
 
    IMPLICIT NONE
@@ -103 +104 @@
       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
       !!----------------------------------------------------------------------
+      IF( nn_timing == 1 )   CALL timing_start('dia_hth')
 
       IF( kt == nit000 ) THEN
@@ -160 +162 @@
          DO ji = 1, jpi
             IF( tmask(ji,jj,nla10) == 1. ) THEN
-               zu  =  1779.50 + 11.250*tn(ji,jj,nla10) - 3.80*sn(ji,jj,nla10) - 0.0745*tn(ji,jj,nla10)*tn(ji,jj,nla10)   &
-                  &                                                           - 0.0100*tn(ji,jj,nla10)*sn(ji,jj,nla10)
-               zv  =  5891.00 + 38.000*tn(ji,jj,nla10) + 3.00*sn(ji,jj,nla10) - 0.3750*tn(ji,jj,nla10)*tn(ji,jj,nla10)
-               zut =    11.25 -  0.149*tn(ji,jj,nla10) - 0.01*sn(ji,jj,nla10)
-               zvt =    38.00 -  0.750*tn(ji,jj,nla10)
+               zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
+                  &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
+                  &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
+               zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
+                  &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
+               zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
+               zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
                zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
                zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
@@ -184 +188 @@
                !
                zzdep = fsdepw(ji,jj,jk)
-               zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
+               zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
                zzdep = zzdep * tmask(ji,jj,1)
 
@@ -221 +225 @@
                zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
                !
-               zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk)                  ! - delta T(10m)
+               zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
                IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
                IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
@@ -254 +258 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zztmp = tn(ji,jj,jk)
+               zztmp = tsn(ji,jj,jk,jp_tem)
                IF( zztmp >= 20. )   ik20(ji,jj) = jk
                IF( zztmp >= 28. )   ik28(ji,jj) = jk
@@ -273 +277 @@
                zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
                   &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
-                  &  * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   &
-                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
+                  &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
+                  &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
                hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
             ELSE 
@@ -284 +288 @@
                zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation
                   &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
-                  &  * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   &
-                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
+                  &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
+                  &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
                hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
             ELSE 
@@ -309 +313 @@
       ! surface boundary condition
       IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
-      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk)   
+      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,jk,jp_tem) * sshn(:,:) * tmask(:,:,jk)   
       ENDIF
       ! integration down to ilevel
       DO jk = 1, ilevel
          zthick(:,:) = zthick(:,:) + fse3t(:,:,jk)
-         htc3  (:,:) = htc3  (:,:) + fse3t(:,:,jk) * tn(:,:,jk) * tmask(:,:,jk)
+         htc3  (:,:) = htc3  (:,:) + fse3t(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
       END DO
       ! deepest layer
@@ -320 +324 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            htc3(ji,jj) = htc3(ji,jj) + tn(ji,jj,ilevel+1) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
+            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )  &
+                                                                   * tmask(ji,jj,ilevel+1)
+            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )   &
+               &                                                   * tmask(ji,jj,ilevel+1)
          END DO
       END DO
@@ -327 +334 @@
       htc3(:,:) = zcoef * htc3(:,:)
       CALL iom_put( "hc300", htc3 )      ! first 300m heat content
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_hth')
       !
    END SUBROUTINE dia_hth

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

@@ -29 +29 @@
    USE lib_mpp          ! MPP library
    USE lbclnk           ! lateral boundary condition - processor exchanges
+   USE timing           ! preformance summary
+   USE wrk_nemo         ! working arrays
 
    IMPLICIT NONE
@@ -95 +97 @@
       !!----------------------------------------------------------------------
       INTEGER               ::   dia_ptr_alloc   ! return value
-      INTEGER, DIMENSION(5) ::   ierr
+      INTEGER, DIMENSION(6) ::   ierr
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -121 +123 @@
          &     ndex_h_ind_30(jpj),   ndex_h_ipc_30(jpj), Stat=ierr(5) )
          !
+     ALLOCATE( btm30(jpi,jpj) , STAT=ierr(6)  )
+         !
       dia_ptr_alloc = MAXVAL( ierr )
       IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc )
@@ -209 +213 @@
       !! ** Action  : - p_fval: i-mean poleward flux of pva
       !!----------------------------------------------------------------------
-#if defined key_mpp_mpi
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zwork => wrk_1d_1
-#endif
       !!
       IMPLICIT none
@@ -225 +225 @@
       INTEGER               ::   ijpjjpk
 #endif
+#if defined key_mpp_mpi
+      REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point
+#endif
       !!--------------------------------------------------------------------
       !
 #if defined key_mpp_mpi
-      IF( wrk_in_use(1, 1) ) THEN
-         CALL ctl_stop('ptr_vjk: ERROR - requested workspace array is unavailable')   ;   RETURN
-      END IF
+      ijpjjpk = jpj*jpk
+      CALL wrk_alloc( jpj*jpk, zwork )
 #endif
 
@@ -257 +259 @@
       !
 #if defined key_mpp_mpi
-      ijpjjpk = jpj*jpk
       ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk
       zwork(1:ijpjjpk) = RESHAPE( p_fval, ish )
@@ -265 +266 @@
       !
 #if defined key_mpp_mpi
-      IF( wrk_not_released(1, 1) )   CALL ctl_stop('ptr_vjk: ERROR - failed to release workspace array')
+      CALL wrk_dealloc( jpj*jpk, zwork )
 #endif
       !
@@ -281 +282 @@
       !! ** Action  : - p_fval: i-sum of e1t*e3t*pta
       !!----------------------------------------------------------------------
-#if defined key_mpp_mpi
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zwork => wrk_1d_1
-#endif
       !!
       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! tracer flux array at T-point
@@ -296 +293 @@
       INTEGER               ::   ijpjjpk
 #endif
+#if defined key_mpp_mpi
+      REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point
+#endif
       !!-------------------------------------------------------------------- 
       !
 #if defined key_mpp_mpi
-      IF( wrk_in_use(1, 1) ) THEN
-         CALL ctl_stop('ptr_tjk: requested workspace array unavailable')   ;   RETURN
-      ENDIF
+      ijpjjpk = jpj*jpk
+      CALL wrk_alloc( jpj*jpk, zwork )
 #endif
 
@@ -315 +314 @@
       END DO
 #if defined key_mpp_mpi
-      ijpjjpk = jpj*jpk
       ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk
       zwork(1:ijpjjpk)= RESHAPE( p_fval, ish )
@@ -323 +321 @@
       !
 #if defined key_mpp_mpi
-      IF( wrk_not_released(1, 1) )   CALL ctl_stop('ptr_tjk: failed to release workspace array')
+      CALL wrk_dealloc( jpj*jpk, zwork )
 #endif
       !    
@@ -343 +341 @@
       !!----------------------------------------------------------------------
       !
+      IF( nn_timing == 1 )   CALL timing_start('dia_ptr')
+      !
       IF( kt == nit000 .OR. MOD( kt, nn_fptr ) == 0 )   THEN
          !
@@ -349 +349 @@
             IF( ln_diaznl ) THEN               ! i-mean temperature and salinity
                DO jn = 1, nptr
-                  tn_jk(:,:,jn) = ptr_tjk( tn(:,:,:), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
+                  tn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
                END DO
             ENDIF
@@ -368 +368 @@
             !
             !                          ! Transports
-            !                                ! local heat & salt transports at T-points  ( tn*mj[vn+v_eiv] )
+            !                                ! local heat & salt transports at T-points  ( tsn*mj[vn+v_eiv] )
             vt(:,:,jpk) = 0._wp   ;   vs(:,:,jpk) = 0._wp
             DO jk= 1, jpkm1
@@ -378 +378 @@
                      zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) ) * 0.5_wp
 #endif 
-                     vt(:,jj,jk) = zv * tn(:,jj,jk)
-                     vs(:,jj,jk) = zv * sn(:,jj,jk)
+                     vt(:,jj,jk) = zv * tsn(:,jj,jk,jp_tem)
+                     vs(:,jj,jk) = zv * tsn(:,jj,jk,jp_sal)
                   END DO
                END DO
@@ -430 +430 @@
       ENDIF
       !
-      IF( kt == nitend )   CALL histclo( numptr )      ! Close the file
+#if defined key_mpp_mpi
+      IF( kt == nitend .AND. l_znl_root )   CALL histclo( numptr )      ! Close the file
+#else
+      IF( kt == nitend )                    CALL histclo( numptr )      ! Close the file
+#endif
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr')
       !
    END SUBROUTINE dia_ptr
@@ -449 +455 @@
       NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, ln_ptrcomp, nn_fptr, nn_fwri
       !!----------------------------------------------------------------------
-
-      !                                      ! allocate dia_ptr arrays
-      IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate arrays' )
+      IF( nn_timing == 1 )   CALL timing_start('dia_ptr_init')
 
       REWIND( numnam )                 ! Read Namelist namptr : poleward transport parameters
@@ -468 +472 @@
          WRITE(numout,*) '      Frequency of outputs                               nn_fwri    = ', nn_fwri
       ENDIF
-
-      IF( ln_subbas ) THEN   ;   nptr = 5       ! Global, Atlantic, Pacific, Indian, Indo-Pacific
-      ELSE                   ;   nptr = 1       ! Global only
-      ENDIF
-
-      rc_pwatt = rc_pwatt * rau0 * rcp          ! conversion from K.s-1 to PetaWatt
-
-      IF( .NOT. ln_diaptr ) THEN       ! diaptr not used
-        RETURN
-      ENDIF
       
-      IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
-
-      IF( ln_subbas ) THEN                ! load sub-basin mask
-         CALL iom_open( 'subbasins', inum )
-         CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
-         CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
-         CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
-         CALL iom_close( inum )
-         btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
-         WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
-         ELSE WHERE                     ;   btm30(:,:) = tmask(:,:,1)
-         END WHERE
-      ENDIF
-      btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
+      IF( ln_diaptr) THEN  
       
-      DO jn = 1, nptr
-         btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
-      END DO
+         IF( ln_subbas ) THEN   ;   nptr = 5       ! Global, Atlantic, Pacific, Indian, Indo-Pacific
+         ELSE                   ;   nptr = 1       ! Global only
+         ENDIF
+
+         !                                      ! allocate dia_ptr arrays
+         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
+
+         rc_pwatt = rc_pwatt * rau0 * rcp          ! conversion from K.s-1 to PetaWatt
+
+         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum
+
+         IF( ln_subbas ) THEN                ! load sub-basin mask
+            CALL iom_open( 'subbasins', inum )
+            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin
+            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin
+            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin
+            CALL iom_close( inum )
+            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin
+            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
+            ELSE WHERE                     ;   btm30(:,:) = tmask(:,:,1)
+            END WHERE
+         ENDIF
+         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean
       
-      IF( lk_vvl )   CALL ctl_stop( 'diaptr: error in vvl case as constant i-mean surface is used' )
-
-      !                                   ! i-sum of e1v*e3v surface and its inverse
-      DO jn = 1, nptr
-         sjk(:,:,jn) = ptr_tjk( tmask(:,:,:), btmsk(:,:,jn) )
-         r1_sjk(:,:,jn) = 0._wp
-         WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
-      END DO
+         DO jn = 1, nptr
+            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only
+         END DO
+      
+         IF( lk_vvl )   CALL ctl_stop( 'diaptr: error in vvl case as constant i-mean surface is used' )
+
+         !                                   ! i-sum of e1v*e3v surface and its inverse
+         DO jn = 1, nptr
+            sjk(:,:,jn) = ptr_tjk( tmask(:,:,:), btmsk(:,:,jn) )
+            r1_sjk(:,:,jn) = 0._wp
+            WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn)
+         END DO
+
+      ! Initialise arrays to zero because diatpr is called before they are first calculated
+      ! Note that this means diagnostics will not be exactly correct when model run is restarted.
+      htr_adv(:) = 0._wp ; str_adv(:) =  0._wp ;  htr_ldf(:) = 0._wp ; str_ldf(:) =  0._wp
 
 #if defined key_mpp_mpi 
-      iglo (1) = jpjglo                   ! MPP case using MPI  ('key_mpp_mpi')
-      iloc (1) = nlcj
-      iabsf(1) = njmppt(narea)
-      iabsl(:) = iabsf(:) + iloc(:) - 1
-      ihals(1) = nldj - 1
-      ihale(1) = nlcj - nlej
-      idid (1) = 2
-      CALL flio_dom_set( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_ptr )
+         iglo (1) = jpjglo                   ! MPP case using MPI  ('key_mpp_mpi')
+         iloc (1) = nlcj
+         iabsf(1) = njmppt(narea)
+         iabsl(:) = iabsf(:) + iloc(:) - 1
+         ihals(1) = nldj - 1
+         ihale(1) = nlcj - nlej
+         idid (1) = 2
+         CALL flio_dom_set( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_ptr )
 #else
-      nidom_ptr = FLIO_DOM_NONE
-#endif
+         nidom_ptr = FLIO_DOM_NONE
+#endif
+      ENDIF 
+      ! 
+      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr_init')
       ! 
    END SUBROUTINE dia_ptr_init
@@ -531 +543 @@
       !! ** Method  :   NetCDF file
       !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zphi => wrk_1d_1, zfoo => wrk_1d_2    ! 1D workspace
-      USE wrk_nemo, ONLY:   z_1  => wrk_2d_1                      ! 2D      -
       !!
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
@@ -548 +557 @@
 #endif
       REAL(wp)           ::   zsto, zout, zdt, zjulian                ! temporary scalars
-      !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(1, 1,2) .OR. wrk_in_use(2, 1) ) THEN
-         CALL ctl_stop('dia_ptr_wri: requested workspace arrays unavailable')   ;   RETURN
-      ENDIF
+      !!
+      REAL(wp), POINTER, DIMENSION(:)   ::   zphi, zfoo    ! 1D workspace
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z_1           ! 2D workspace
+      !!-------------------------------------------------------------------- 
+      !
+      CALL wrk_alloc( jpi      , zphi , zfoo )
+      CALL wrk_alloc( jpi , jpk, z_1 )
 
       ! define time axis
@@ -866 +877 @@
       ENDIF
       !
-      IF( wrk_not_released(1, 1,2) .OR.    &
-          wrk_not_released(2, 1)    )   CALL ctl_stop('dia_ptr_wri: failed to release workspace arrays')
+      CALL wrk_dealloc( jpi      , zphi , zfoo )
+      CALL wrk_dealloc( jpi , jpk, z_1 )
       !
   END SUBROUTINE dia_ptr_wri

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

@@ -28 +28 @@
    USE ldftra_oce      ! ocean active tracers: lateral physics
    USE ldfdyn_oce      ! ocean dynamics: lateral physics
+   USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
    USE sol_oce         ! solver variables
    USE sbc_oce         ! Surface boundary condition: ocean fields
@@ -46 +47 @@
    USE limwri_2 
 #endif
-   USE dtatem
-   USE dtasal
    USE lib_mpp         ! MPP library
+   USE timing          ! preformance summary
+   USE wrk_nemo        ! working array
 
    IMPLICIT NONE
@@ -116 +117 @@
       !! ** Method  :  use iom_put
       !!----------------------------------------------------------------------
-      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
-      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY: z2d => wrk_2d_1
       !!
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
@@ -124 +122 @@
       INTEGER                      ::   ji, jj, jk              ! dummy loop indices
       REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
+      !!
+      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d       ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
-      IF( wrk_in_use(2, 1))THEN
-         CALL ctl_stop('dia_wri: ERROR - requested 2D workspace unavailable.')
-         RETURN
-      END IF
+      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
+      ! 
+      CALL wrk_alloc( jpi , jpj      , z2d )
+      CALL wrk_alloc( jpi , jpj, jpk , z3d )
       !
       ! Output the initial state and forcings
@@ -137 +138 @@
       ENDIF
 
-      CALL iom_put( "toce"   , tn                    )    ! temperature
-      CALL iom_put( "soce"   , sn                    )    ! salinity
-      CALL iom_put( "sst"    , tn(:,:,1)             )    ! sea surface temperature
-      CALL iom_put( "sst2"   , tn(:,:,1) * tn(:,:,1) )    ! square of sea surface temperature
-      CALL iom_put( "sss"    , sn(:,:,1)             )    ! sea surface salinity
-      CALL iom_put( "sss2"   , sn(:,:,1) * sn(:,:,1) )    ! square of sea surface salinity
-      CALL iom_put( "uoce"   , un                    )    ! i-current      
-      CALL iom_put( "voce"   , vn                    )    ! j-current
+      CALL iom_put( "toce"   , tsn(:,:,:,jp_tem)                     )    ! temperature
+      CALL iom_put( "soce"   , tsn(:,:,:,jp_sal)                     )    ! salinity
+      CALL iom_put( "sst"    , tsn(:,:,1,jp_tem)                     )    ! sea surface temperature
+      CALL iom_put( "sst2"   , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) )    ! square of sea surface temperature
+      CALL iom_put( "sss"    , tsn(:,:,1,jp_sal)                     )    ! sea surface salinity
+      CALL iom_put( "sss2"   , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) )    ! square of sea surface salinity
+      CALL iom_put( "uoce"   , un                                    )    ! i-current      
+      CALL iom_put( "voce"   , vn                                    )    ! j-current
       
-      CALL iom_put( "avt"    , avt                   )    ! T vert. eddy diff. coef.
-      CALL iom_put( "avm"    , avmu                  )    ! T vert. eddy visc. coef.
+      CALL iom_put( "avt"    , avt                                   )    ! T vert. eddy diff. coef.
+      CALL iom_put( "avm"    , avmu                                  )    ! T vert. eddy visc. coef.
       IF( lk_zdfddm ) THEN
-         CALL iom_put( "avs" , fsavs(:,:,:)          )    ! S vert. eddy diff. coef.
+         CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef.
       ENDIF
 
       DO jj = 2, jpjm1                                    ! sst gradient
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            zztmp      = tn(ji,jj,1)
-            zztmpx     = ( tn(ji+1,jj  ,1) - zztmp ) / e1u(ji,jj) + ( zztmp - tn(ji-1,jj  ,1) ) / e1u(ji-1,jj  )
-            zztmpy     = ( tn(ji  ,jj+1,1) - zztmp ) / e2v(ji,jj) + ( zztmp - tn(ji  ,jj-1,1) ) / e2v(ji  ,jj-1)
+            zztmp      = tsn(ji,jj,1,jp_tem)
+            zztmpx     = ( tsn(ji+1,jj  ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  )
+            zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / e2v(ji  ,jj-1)
             z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
                &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
@@ -178 +179 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
                END DO
             END DO
@@ -192 +193 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) )
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
                END DO
             END DO
@@ -200 +201 @@
       ENDIF
       !
-      IF( wrk_not_released(2, 1))THEN
-         CALL ctl_stop('dia_wri: ERROR - failed to release 2D workspace.')
-         RETURN
-      END IF
+      CALL wrk_dealloc( jpi , jpj      , z2d )
+      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
       !
    END SUBROUTINE dia_wri
@@ -224 +225 @@
       !!      Each nwrite time step, output the instantaneous or mean fields
       !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY: zw2d => wrk_2d_1
       !!
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
@@ -232 +231 @@
       CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
       INTEGER  ::   inum = 11                                ! temporary logical unit
+      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
+      INTEGER  ::   ierr                                     ! error code return from allocation
       INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
       REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
+      !!
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
       !!----------------------------------------------------------------------
-      !
-      IF( wrk_in_use(2, 1))THEN
-         CALL ctl_stop('dia_wri: ERROR - requested 2D workspace unavailable.')
-         RETURN
-      END IF
+      ! 
+      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
+      !
+      CALL wrk_alloc( jpi , jpj      , zw2d )
+      IF ( ln_traldf_gdia )  call wrk_alloc( jpi , jpj , jpk  , zw3d )
       !
       ! Output the initial state and forcings
@@ -446 +450 @@
          CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
             &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
+         IF( ln_traldf_gdia ) THEN
+            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
+                 &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
+         ELSE
 #if defined key_diaeiv
-         CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
+            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
             &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
 #endif
+         END IF
          !                                                                                      !!! nid_U : 2D
          CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
@@ -459 +468 @@
          CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
             &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
+         IF( ln_traldf_gdia ) THEN
+            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
+                 &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
+         ELSE 
 #if defined key_diaeiv
-         CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
+            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
             &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
 #endif
+         END IF
          !                                                                                      !!! nid_V : 2D
          CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
@@ -472 +486 @@
          CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
             &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
+         IF( ln_traldf_gdia ) THEN
+            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
+                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
+         ELSE
 #if defined key_diaeiv
-         CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
-            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
-#endif
+            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
+                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
+#endif
+         END IF
          CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
             &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
@@ -516 +535 @@
 
       ! Write fields on T grid
-      CALL histwrite( nid_T, "votemper", it, tn            , ndim_T , ndex_T  )   ! temperature
-      CALL histwrite( nid_T, "vosaline", it, sn            , ndim_T , ndex_T  )   ! salinity
-      CALL histwrite( nid_T, "sosstsst", it, tn(:,:,1)     , ndim_hT, ndex_hT )   ! sea surface temperature
-      CALL histwrite( nid_T, "sosaline", it, sn(:,:,1)     , ndim_hT, ndex_hT )   ! sea surface salinity
+      CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T  )   ! temperature
+      CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T  )   ! salinity
+      CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem), ndim_hT, ndex_hT )   ! sea surface temperature
+      CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal), ndim_hT, ndex_hT )   ! sea surface salinity
       CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
 !!$#if  defined key_lim3 || defined key_lim2 
@@ -528 +547 @@
 !!$      CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff
       CALL histwrite( nid_T, "sowaflcd", it, ( emps-rnf )  , ndim_hT, ndex_hT )   ! c/d water flux
-      zw2d(:,:) = ( emps(:,:) - rnf(:,:) ) * sn(:,:,1) * tmask(:,:,1)
+      zw2d(:,:) = ( emps(:,:) - rnf(:,:) ) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
       CALL histwrite( nid_T, "sosalflx", it, zw2d          , ndim_hT, ndex_hT )   ! c/d salt flux
       CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
@@ -539 +558 @@
       CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
       CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-      IF( ln_ssr ) zw2d(:,:) = erp(:,:) * sn(:,:,1) * tmask(:,:,1)
+      IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
       CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
 #endif
@@ -545 +564 @@
       CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
       CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * sn(:,:,1) * tmask(:,:,1)
+         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
       CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
 #endif
@@ -570 +589 @@
          ! Write fields on U grid
       CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
+      IF( ln_traldf_gdia ) THEN
+         IF (.not. ALLOCATED(psix_eiv))THEN
+            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
+            IF( lk_mpp   )   CALL mpp_sum ( ierr )
+            IF( ierr > 0 )   CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
+            psix_eiv(:,:,:) = 0.0_wp
+            psiy_eiv(:,:,:) = 0.0_wp
+         ENDIF
+         DO jk=1,jpkm1
+            zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
+         END DO
+         zw3d(:,:,jpk) = 0._wp
+         CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U )           ! i-eiv current
+      ELSE
 #if defined key_diaeiv
-      CALL histwrite( nid_U, "vozoeivu", it, u_eiv         , ndim_U , ndex_U )    ! i-eiv current
-#endif
+         CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U )          ! i-eiv current
+#endif
+      ENDIF
       CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
 
          ! Write fields on V grid
       CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
+      IF( ln_traldf_gdia ) THEN
+         DO jk=1,jpk-1
+            zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
+         END DO
+         zw3d(:,:,jpk) = 0._wp
+         CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V )           ! j-eiv current
+      ELSE
 #if defined key_diaeiv
-      CALL histwrite( nid_V, "vomeeivv", it, v_eiv         , ndim_V , ndex_V  )   ! j-eiv current
-#endif
+         CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V )          ! j-eiv current
+#endif
+      ENDIF
       CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
          ! Write fields on W grid
       CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
+      IF( ln_traldf_gdia ) THEN
+         DO jk=1,jpk-1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1  ! vector opt.
+                  zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
+                       &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
+               END DO
+            END DO
+         END DO
+         zw3d(:,:,jpk) = 0._wp
+         CALL histwrite( nid_W, "voveeivw", it, zw3d          , ndim_T, ndex_T )    ! vert. eiv current
+      ELSE
 #   if defined key_diaeiv
-      CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
+         CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
 #   endif
+      ENDIF
       CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
       CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
@@ -608 +663 @@
       ENDIF
       !
-      IF( wrk_not_released(2, 1))THEN
-         CALL ctl_stop('dia_wri: ERROR - failed to release 2D workspace.')
-         RETURN
-      END IF
+      CALL wrk_dealloc( jpi , jpj      , zw2d )
+      IF ( ln_traldf_gdia )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
       !
    END SUBROUTINE dia_wri
@@ -640 +695 @@
       REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
       !!----------------------------------------------------------------------
+      ! 
+      IF( nn_timing == 1 )   CALL timing_start('dia_wri_state')
 
       ! 0. Initialisation
@@ -711 +768 @@
 
       ! Write all fields on T grid
-      CALL histwrite( id_i, "votemper", kt, tn      , jpi*jpj*jpk, idex )    ! now temperature
-      CALL histwrite( id_i, "vosaline", kt, sn      , jpi*jpj*jpk, idex )    ! now salinity
-      CALL histwrite( id_i, "sossheig", kt, sshn     , jpi*jpj    , idex )    ! sea surface height
-      CALL histwrite( id_i, "vozocrtx", kt, un       , jpi*jpj*jpk, idex )    ! now i-velocity
-      CALL histwrite( id_i, "vomecrty", kt, vn       , jpi*jpj*jpk, idex )    ! now j-velocity
-      CALL histwrite( id_i, "vovecrtz", kt, wn       , jpi*jpj*jpk, idex )    ! now k-velocity
-      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf), jpi*jpj    , idex )    ! freshwater budget
-      CALL histwrite( id_i, "sohefldo", kt, qsr + qns, jpi*jpj    , idex )    ! total heat flux
-      CALL histwrite( id_i, "soshfldo", kt, qsr      , jpi*jpj    , idex )    ! solar heat flux
-      CALL histwrite( id_i, "soicecov", kt, fr_i     , jpi*jpj    , idex )    ! ice fraction
-      CALL histwrite( id_i, "sozotaux", kt, utau     , jpi*jpj    , idex )    ! i-wind stress
-      CALL histwrite( id_i, "sometauy", kt, vtau     , jpi*jpj    , idex )    ! j-wind stress
+      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
+      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
+      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
+      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
+      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
+      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
+      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf )       , jpi*jpj    , idex )    ! freshwater budget
+      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
+      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
+      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
+      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
+      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
 
       ! 3. Close the file
@@ -735 +792 @@
       ENDIF
 #endif
+       
+      IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state')
+      ! 
 
    END SUBROUTINE dia_wri_state

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri_dimg.h90

@@ -82 +82 @@
     INTEGER :: inbsel, jk
     INTEGER :: iyear,imon,iday
+    INTEGER :: ialloc
     REAL(wp) :: zdtj
     CHARACTER(LEN=80) :: clname
@@ -88 +89 @@
     CHARACTER(LEN= 4) :: clver
     !!----------------------------------------------------------------------
+    IF( nn_timing == 1 )   CALL timing_start('dia_wri')
     !
     !  Initialization
     !  ---------------
     !
-    IF(.not.ALLOCATED(um))THEN
+    IF( .NOT. ALLOCATED(um) )THEN
        ALLOCATE(um(jpi,jpj,jpk), vm(jpi,jpj,jpk), &
                 wm(jpi,jpj,jpk),                  &
@@ -98 +100 @@
                 tm(jpi,jpj,jpk), sm(jpi,jpj,jpk), &
                 fsel(jpi,jpj,jpk),                &
-                Stat=jk)
-       IF(jk /= 0)THEN
-          WRITE(*,*) 'ERROR: allocate failed in dia_wri (diawri_dimg.h90)'
-          CALL mppabort()
-       END IF
-    END IF
+                STAT=ialloc )
+       !
+       IF( lk_mpp      )   CALL mpp_sum ( ialloc  )
+       IF( ialloc /= 0 )   CALL ctl_warn('dia_wri( diawri_dimg.h90) : failed to allocate arrays')
+    ENDIF
+
 
     inbsel = 17
@@ -152 +154 @@
        wm(:,:,:)=wm(:,:,:) + wn (:,:,:)
        avtm(:,:,:)=avtm(:,:,:) + avt (:,:,:)
-       tm(:,:,:)=tm(:,:,:) + tn (:,:,:)
-       sm(:,:,:)=sm(:,:,:) + sn (:,:,:)
+       tm(:,:,:)=tm(:,:,:) + tsn(:,:,:,jp_tem)
+       sm(:,:,:)=sm(:,:,:) + tsn(:,:,:,jp_sal)
        !
        fsel(:,:,1 ) = fsel(:,:,1 ) + utau(:,:) * umask(:,:,1)
@@ -159 +161 @@
        fsel(:,:,3 ) = fsel(:,:,3 ) + qsr (:,:) + qns  (:,:) 
        fsel(:,:,4 ) = fsel(:,:,4 ) + ( emp(:,:)-rnf(:,:) ) 
-       !        fsel(:,:,5 ) = fsel(:,:,5 ) + tb  (:,:,1)  !RB not used
+       !        fsel(:,:,5 ) = fsel(:,:,5 ) + tsb(:,:,1,jp_tem)  !RB not used
        fsel(:,:,6 ) = fsel(:,:,6 ) + sshn(:,:) 
        fsel(:,:,7 ) = fsel(:,:,7 ) + qsr(:,:)
@@ -226 +228 @@
           fsel(:,:,3 ) = (qsr (:,:) + qns (:,:)) * tmask(:,:,1)
           fsel(:,:,4 ) = ( emp(:,:)-rnf(:,:) ) * tmask(:,:,1) 
-          !         fsel(:,:,5 ) = (tb  (:,:,1) - sf_sst(1)%fnow(:,:) ) *tmask(:,:,1) !RB not used
+          !         fsel(:,:,5 ) = (tsb(:,:,1,jp_tem) - sf_sst(1)%fnow(:,:) ) *tmask(:,:,1) !RB not used
 
           fsel(:,:,6 ) = sshn(:,:)
@@ -302 +304 @@
 
        IF( ll_dia_inst) THEN
-          CALL dia_wri_dimg(clname, cltext, tn, jpk, 'T')
-       ELSE
-          CALL dia_wri_dimg(clname, cltext, tm, jpk, 'T')
+          CALL dia_wri_dimg(clname, cltext, tsn(:,:,:,jp_tem), jpk, 'T')
+       ELSE
+          CALL dia_wri_dimg(clname, cltext, tm               , jpk, 'T')
        ENDIF
        !
@@ -314 +316 @@
 
        IF( ll_dia_inst) THEN
-          CALL dia_wri_dimg(clname, cltext, sn, jpk, 'T')
-       ELSE
-          CALL dia_wri_dimg(clname, cltext, sm, jpk, 'T')
+          CALL dia_wri_dimg(clname, cltext, tsn(:,:,:,jp_sal), jpk, 'T')
+       ELSE
+          CALL dia_wri_dimg(clname, cltext, sm               , jpk, 'T')
        ENDIF
        !
@@ -357 +359 @@
     ENDIF
     !
+    IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
+    !
 9000 FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc")
     !