source: trunk/NEMO/OPA_SRC/DIA/diawri_dimg.h90 @ 642

  !!----------------------------------------------------------------------
  !!                        ***  diawri_dimg.h90  ***
  !!----------------------------------------------------------------------
  !!   OPA 9.0 , LOCEAN-IPSL (2005) 
  !! $Header$ 
  !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
  !!----------------------------------------------------------------------

  SUBROUTINE dia_wri (kt, kindic)
    !!----------------------------------------------------------------------
    !!           *** routine dia_wri ***
    !!
    !! ** Purpose : output dynamics and tracer fields on direct access file
    !!              suitable for MPP computing
    !!
    !! ** define key : 'key_dimgout'
    !!
    !! **  Method : Default is to cumulate the values over the interval between
    !!      2 output, and each nwrite time-steps  the mean value is  computed 
    !!      and written to the direct access file.
    !!     If 'key_diainstant' is defined, no mean values are computed and the 
    !!     instantaneous fields are dump.
    !!     kindic is 0 or >0 in normal condition. When < 0 it indicates an error
    !!     condition and instantaneous file output is forced.
    !!       Each processor creates its own file with its local data
    !!     Merging all the files is performed off line by a dedicated program
    !!
    !! ** Arguments :
    !!     kt      : time-step number
    !!     kindinc :  error condition indicator : >=0 :  OK, < 0 : error.
    !!
    !! ** Naming convention for files
    !!
    !! {cexper}_{var}_y----m--d--.dimg
    !!   cexper is the name of the experience, given in the namelist
    !!   var can be either U, V, T, S, KZ, SSH, ...
    !!   var can also be 2D, which means that each level of the file is a 2D field as described below
    !!    y----m--d--  is the date at the time of the dump
    !!    For mpp output, each processor dumps its own memory, on appropriate record range 
    !!    (direct access : for level jk of a klev field on proc narea irec = 1+ klev*(narea -1) + jk )
    !!    To be tested with a lot of procs !!!!
    !!
    !!  level 1:  taux(:,:) * umask(:,:,1) zonal stress in N.m-2
    !!  level 2:  tauy(:,:) * vmask(:,:,1) meridional stress in N. m-2
    !!  level 3:   qt  (:,:)               total heat flux (W/m2)
    !!  level 4:   emp (:,:)               E-P flux (mm/day)
    !!  level 5:  tb  (:,:,1)-sst          model SST -forcing sst (degree C)
    !!  level 6:  bsfb(:,:)         streamfunction (m**3/s)
    !!  level 7:  qsr (:,:)         solar flux (W/m2)
    !!  level 8:  qrp (:,:)                relax component of T flux.
    !!  level 9:  erp (:,:)                relax component of S flux
    !!  level 10: hmld(:,:)                turbocline depth
    !!  level 11: hmlp(:,:)                mixed layer depth
    !!  level 12: freeze (:,:)             Ice cover (1. or 0.)
    !!  level 13: sst(:,:)                 the observed SST we relax to.
    !!  level 14: qct(:,:)                 equivalent flux due to treshold SST
    !!  level 15: fbt(:,:)                 feedback term .
    !!  level 16: emps(:,:)                concentration/dilution water flux
    !!  level 17: fsalt(:,:)               Ice=>ocean net freshwater
    !!  level 18: gps(:,:)                 the surface pressure (m).
    !!  level 19: spgu(:,:)                the surface pressure gradient in X direction.
    !!  level 20: spgv(:,:)                the surface pressure gradient in Y direction.
    !! 
    !! History 
    !!      original  : 91-03 ()
    !!      additions : 91-11 (G. Madec)
    !!      additions : 92-06 (M. Imbard) correction restart file
    !!      additions : 92-07 (M. Imbard) split into diawri and rstwri
    !!      additions : 93-03 (M. Imbard) suppress writibm
    !!      additions : 94-12 (M. Imbard) acces direct files
    !!      additions : 97-2002 ( Clipper Group ) dimg files
    !!                  dec 2003 ( J.M. Molines) f90, mpp output for OPA9.0
    !!   9.0  !  05-05  (S. Theetten) add emps fsalt move gps spgu spgv 2 lines below
    !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
    !!----------------------------------------------------------------------
    !! * modules used
    USE lib_mpp
    USE dtasst, ONLY : sst

    !! * Arguments
    INTEGER ,INTENT(in) :: kt, kindic

    !! * local declarations
    INTEGER :: inbsel, jk
!!  INTEGER :: iwrite
    INTEGER :: iyear,imon,iday
    INTEGER, SAVE :: nmoyct 

#if defined key_diainstant
    LOGICAL, PARAMETER :: ll_dia_inst=.TRUE.  !: for instantaneous output
#else
    LOGICAL, PARAMETER :: ll_dia_inst=.FALSE. !: for average output
#endif

    REAL(wp), SAVE, DIMENSION (jpi,jpj,jpk) ::  um , vm   ! used to compute mean u, v fields
    REAL(wp), SAVE, DIMENSION (jpi,jpj,jpk) ::  wm        ! used to compute mean w fields
    REAL(wp), SAVE, DIMENSION (jpi,jpj,jpk) ::  avtm      ! used to compute mean kz fields
    REAL(wp), SAVE, DIMENSION (jpi,jpj,jpk) ::  tm , sm   ! used to compute mean t, s fields
    REAL(wp), SAVE, DIMENSION (jpi,jpj,jpk) ::  fsel      ! used to compute mean 2d fields
    REAL(wp) :: zdtj
    !
    CHARACTER(LEN=80) :: clname
    CHARACTER(LEN=80) :: cltext
    CHARACTER(LEN=80) :: clmode
    CHARACTER(LEN= 4) :: clver
    !
    !  Initialization
    !  ---------------
    !
#ifdef key_diaspr
    inbsel = 20
#else
    inbsel = 17
#endif
#if defined key_flx_core
    inbsel = 23
#endif

    IF( inbsel >  jpk) THEN
       IF( lwp) WRITE(numout,*)  &
            ' STOP inbsel =',inbsel,' is larger than jpk=',jpk
       STOP
    ENDIF


    iyear = ndastp/10000
    imon = (ndastp-iyear*10000)/100
    iday = ndastp - imon*100 - iyear*10000

    !     
    !! dimg format V1.0 should start with the 4 char. string '@!01'
    !!
    clver='@!01'
    !
    IF( .NOT. ll_dia_inst ) THEN
       !
       !! * Mean output section
       !! ----------------------
       !
       IF( kt == nit000 .AND. lwp ) WRITE(numout,*) &
            'THE OUTPUT FILES CONTAINS THE AVERAGE OF EACH FIELD'
       !
       IF( kt == nit000  ) THEN
          ! reset arrays for average computation
          nmoyct = 0
          !
          um(:,:,:) = 0._wp
          vm(:,:,:) = 0._wp
          wm(:,:,:) = 0._wp
          avtm(:,:,:) = 0._wp
          tm(:,:,:) = 0._wp
          sm(:,:,:) = 0._wp
          fsel(:,:,:) = 0._wp
          !
       ENDIF

       !  cumulate values
       !  ---------------

       nmoyct = nmoyct+1
       !
       um(:,:,:)=um(:,:,:) + un (:,:,:)
       vm(:,:,:)=vm(:,:,:) + vn (:,:,:)
       wm(:,:,:)=wm(:,:,:) + wn (:,:,:)
       avtm(:,:,:)=avtm(:,:,:) + avt (:,:,:)
       tm(:,:,:)=tm(:,:,:) + tn (:,:,:)
       sm(:,:,:)=sm(:,:,:) + sn (:,:,:)
       !
       fsel(:,:,1 ) = fsel(:,:,1 ) + taux(:,:) * umask(:,:,1)
       fsel(:,:,2 ) = fsel(:,:,2 ) + tauy(:,:) * vmask(:,:,1)
       fsel(:,:,3 ) = fsel(:,:,3 ) + qt  (:,:) 
       fsel(:,:,4 ) = fsel(:,:,4 ) + emp (:,:)
       fsel(:,:,5 ) = fsel(:,:,5 ) + tb  (:,:,1) - sst(:,:)
#if ! defined key_dynspg_rl
       fsel(:,:,6 ) = fsel(:,:,6 ) + sshn(:,:)    ! SSH
#else
       fsel(:,:,6 ) = fsel(:,:,6 ) + bsfn(:,:)    ! BSF
#endif
       fsel(:,:,7 ) = fsel(:,:,7 ) + qsr(:,:)
       fsel(:,:,8 ) = fsel(:,:,8 ) + qrp (:,:)
       fsel(:,:,9 ) = fsel(:,:,9 ) + erp (:,:)
       fsel(:,:,10) = fsel(:,:,10) + hmld(:,:)
       fsel(:,:,11) = fsel(:,:,11) + hmlp(:,:)
       fsel(:,:,12) = fsel(:,:,12) + freeze(:,:)
       fsel(:,:,13) = fsel(:,:,13) + sst(:,:)
       !        fsel(:,:,14) = fsel(:,:,14) + qct(:,:)
       !        fsel(:,:,15) = fsel(:,:,15) + fbt(:,:)
       fsel(:,:,16) = fsel(:,:,16) + emps(:,:)
#if defined key_ice_lim
       fsel(:,:,17) = fsel(:,:,17) + fsalt(:,:)
#endif
#ifdef key_diaspr   
       fsel(:,:,18) = fsel(:,:,18) + gps(:,:)/g 
#endif
#if defined key_flx_core
       fsel(:,:,21) = fsel(:,:,21) + qla(:,:)
       fsel(:,:,22) = fsel(:,:,22) + qlw(:,:)
       fsel(:,:,23) = fsel(:,:,23) + qsb(:,:)
#endif
       !
       ! Output of dynamics and tracer fields and selected fields (numwri)
       ! -----------------------------------------------------------------
       !
       !
       zdtj=rdt/86400.   ! time step in days
       WRITE(clmode,'(f5.1,a)' ) nwrite*zdtj,' days average'

       !       iwrite=NINT(adatrj/rwrite)
       !      IF (abs(adatrj-iwrite*rwrite) < zdtj/2.      &

       IF(  ( MOD (kt-nit000+1,nwrite) ==  0 )          &
            &   .OR.       kindic <   0            &
            &   .OR. ( kt == 1 .AND. kindic > 0)  ) THEN
          ! it is time to make a dump on file 
          ! compute average
          um(:,:,:) = um(:,:,:) / nmoyct
          vm(:,:,:) = vm(:,:,:) / nmoyct
          wm(:,:,:) = wm(:,:,:) / nmoyct
          avtm(:,:,:) = avtm(:,:,:) / nmoyct
          tm(:,:,:) = tm(:,:,:) / nmoyct
          sm(:,:,:) = sm(:,:,:) / nmoyct
          !
          fsel(:,:,:) = fsel(:,:,:) / nmoyct
          !
          ! note : the surface pressure is not averaged, but rather 
          ! computed from the averaged gradients.
          !
#ifdef key_diaspr
          fsel(:,:,18)= gps(:,:)/g
          fsel(:,:,19)= spgu(:,:)
          fsel(:,:,20)= spgv(:,:)
#endif
          ! mask mean field with tmask except taux tauy (1,2)
          DO jk=3,inbsel
            fsel(:,:,jk)=fsel(:,:,jk)*tmask(:,:,1)
          END DO
       ENDIF
       !
    ELSE   ! ll_dia_inst true
       !
       !! * Instantaneous output section
       !! ------------------------------
       !
       IF( kt == nit000 .AND. lwp ) WRITE(numout,*) &
            'THE OUTPUT FILES CONTAINS INSTANTANEOUS VALUES OF EACH FIELD'
       !
       zdtj=rdt/86400.   ! time step in days
       !  iwrite=NINT(adatrj/rwrite)
       clmode='instantaneous'
       !     IF (abs(adatrj-iwrite*rwrite) < zdtj/2.  &
       IF (  ( MOD (kt-nit000+1,nwrite) ==  0 )          &
            &   .OR.       kindic <   0            &
            &   .OR. ( kt == 1 .AND. kindic > 0)  ) THEN
          !
          ! transfer wp arrays to sp arrays for dimg files
          fsel(:,:,:) = 0._wp
          !
          fsel(:,:,1 ) = taux(:,:) * umask(:,:,1)
          fsel(:,:,2 ) = tauy(:,:) * vmask(:,:,1)
          fsel(:,:,3 ) = qt  (:,:) * tmask(:,:,1)
          fsel(:,:,4 ) = emp (:,:) * tmask(:,:,1)
          fsel(:,:,5 ) = (tb  (:,:,1) -sst(:,:)) *tmask(:,:,1)

#if ! defined key_dynspg_rl
          fsel(:,:,6 ) = sshn(:,:)
#else
          fsel(:,:,6 ) = bsfn(:,:)
#endif
          fsel(:,:,7 ) = qsr (:,:) * tmask(:,:,1)
          fsel(:,:,8 ) = qrp (:,:) * tmask(:,:,1)
          fsel(:,:,9 ) = erp (:,:) * tmask(:,:,1)
          fsel(:,:,10) = hmld(:,:) * tmask(:,:,1)
          fsel(:,:,11) = hmlp(:,:) * tmask(:,:,1)
          fsel(:,:,12) = freeze(:,:) * tmask(:,:,1)
          fsel(:,:,13) =  sst(:,:)  
          !         fsel(:,:,14) =  qct(:,:)
          !         fsel(:,:,15) =  fbt(:,:)
          fsel(:,:,16) =  emps(:,:) * tmask(:,:,1)
#if defined key_ice_lim
          fsel(:,:,17) =  fsalt(:,:) * tmask(:,:,1)
#endif
#ifdef key_diaspr           
          fsel(:,:,18) =      gps(:,:) /g
          fsel(:,:,19) =      spgu(:,:)
          fsel(:,:,20) =      spgv(:,:)
#endif
#if defined key_flx_core
          fsel(:,:,21) =  qla(:,:)* tmask(:,:,1)
          fsel(:,:,22) =  qlw(:,:)* tmask(:,:,1)
          fsel(:,:,23) =  qsb(:,:)* tmask(:,:,1)
#endif
          !
          !         qct(:,:) = 0._wp
       ENDIF
    ENDIF
    !
    ! Opening of the datrj.out file with the absolute time in day of each dump
    ! this file gives a record of the dump date for post processing ( ASCII file )
    !
    IF(  ( MOD (kt-nit000+1,nwrite) ==  0 )          &
         &   .OR.       kindic <   0            &
         &   .OR. ( kt == 1 .AND. kindic > 0)  ) THEN
       OPEN (14,FILE='datrj.out',FORM='FORMATTED', STATUS='UNKNOWN',POSITION='APPEND ')

       IF( lwp) WRITE(14,'(f10.4,1x,i8)') adatrj, ndastp
       CLOSE(14)

       IF( lwp) WRITE(numout,*)'Days since the begining of the run :',adatrj

       !! * U section

       WRITE(clname,9000) TRIM(cexper),'U',iyear,imon,iday
       cltext=TRIM(cexper)//' U(m/s) '//TRIM(clmode)
       IF ( kindic < 0 )   cltext=TRIM(cexper)//' U(m/s)  instantaneous (explosion)'
       !
       IF( ll_dia_inst) THEN 
          CALL dia_wri_dimg(clname, cltext, un, jpk, 'T')

       ELSE 
          IF( kindic ==  -3 ) THEN
             ! ... in case of explosion on umax, dump instantateous u field instead of mean.
             CALL dia_wri_dimg(clname, cltext, un, jpk, 'T')
          ELSE
             CALL dia_wri_dimg(clname, cltext, um, jpk, 'T')
          ENDIF
       ENDIF

       !! * V section

       WRITE(clname,9000) TRIM(cexper),'V',iyear,imon,iday
       cltext=TRIM(cexper)//' V(m/s) '//TRIM(clmode)
       !
       IF( ll_dia_inst) THEN
          CALL dia_wri_dimg(clname, cltext, vn, jpk, 'T')
       ELSE
          CALL dia_wri_dimg(clname, cltext, vm, jpk, 'T')
       ENDIF
       !

       !! * KZ section

       WRITE(clname,9000) TRIM(cexper),'KZ',iyear,imon,iday
       cltext=TRIM(cexper)//' KZ(m2/s) '//TRIM(clmode)

       IF( ll_dia_inst) THEN
          CALL dia_wri_dimg(clname, cltext, avt, jpk, 'W')
       ELSE
          CALL dia_wri_dimg(clname, cltext, avtm, jpk, 'W')
       ENDIF
       !

       !! * W section

       WRITE(clname,9000) TRIM(cexper),'W',iyear,imon,iday
       cltext=TRIM(cexper)//' W(m/s) '//TRIM(clmode)

       IF( ll_dia_inst) THEN
          CALL dia_wri_dimg(clname, cltext, wn, jpk, 'W')
       ELSE
          CALL dia_wri_dimg(clname, cltext, wm, jpk, 'W')
       ENDIF

       !! * T section

       WRITE(clname,9000) TRIM(cexper),'T',iyear,imon,iday
       cltext=TRIM(cexper)//' T (DegC) '//TRIM(clmode)

       IF( ll_dia_inst) THEN
          CALL dia_wri_dimg(clname, cltext, tn, jpk, 'T')
       ELSE
          CALL dia_wri_dimg(clname, cltext, tm, jpk, 'T')
       ENDIF
       !

       !! * S section

       WRITE(clname,9000) TRIM(cexper),'S',iyear,imon,iday
       cltext=TRIM(cexper)//' S (PSU) '//TRIM(clmode)

       IF( ll_dia_inst) THEN
          CALL dia_wri_dimg(clname, cltext, sn, jpk, 'T')
       ELSE
          CALL dia_wri_dimg(clname, cltext, sm, jpk, 'T')
       ENDIF
       !

       !! * 2D section

       WRITE(clname,9000) TRIM(cexper),'2D',iyear,imon,iday
       cltext='2D fields '//TRIM(clmode)

       IF( ll_dia_inst) THEN
          CALL dia_wri_dimg(clname, cltext, fsel, inbsel, '2')
       ELSE
          CALL dia_wri_dimg(clname, cltext, fsel, inbsel, '2')
       ENDIF

       IF( lk_mpp )   CALL mppsync   ! synchronization in mpp

       !! * Log message in numout 

       IF( lwp)WRITE(numout,*) ' '
       IF( lwp)WRITE(numout,*) ' **** WRITE in numwri ',kt

       IF( lwp .AND.        ll_dia_inst) WRITE(numout,*) '    instantaneous fields'
       IF( lwp .AND. .NOT.  ll_dia_inst) WRITE(numout,*) '    average fields with ',nmoyct,'pdt'
       !
       !
       !! * Reset cumulating arrays  and counter to 0 after writing
       !
       IF( .NOT. ll_dia_inst ) THEN
          nmoyct = 0
          !
          um(:,:,:) = 0._wp
          vm(:,:,:) = 0._wp
          wm(:,:,:) = 0._wp
          tm(:,:,:) = 0._wp
          sm(:,:,:) = 0._wp
          fsel(:,:,:) = 0._wp
          avtm(:,:,:) = 0._wp
       ENDIF
    ENDIF
    !
9000 FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc")

  END SUBROUTINE dia_wri

  SUBROUTINE dia_wri_state ( cdfile_name)
    !!-------------------------------------------------------------------
    !!        ***     ROUTINE dia_wri_state  ***
    !!
    !! ** Purpose :   Dummy routine for compatibility with IOIPSL output
    !!
    !! ** History :
    !!      9.O  ! 03-06  (J.M. Molines ) dimgout
    !!--------------------------------------------------------------------
    !! * Arguments
    CHARACTER (len=*), INTENT(in) ::   cdfile_name   ! name of the file created 
    !!--------------------------------------------------------------------

    IF( lwp) WRITE(numout,*) 'dia_wri_state: Dummy call', cdfile_name
    IF( lwp) WRITE(numout,*) '-------------'
    IF( lwp) WRITE(numout,*)

  END SUBROUTINE dia_wri_state