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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9019 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Ignore:
Timestamp:
2017-12-13T15:58:53+01:00 (6 years ago)
Author:
timgraham
Message:

Merge of dev_CNRS_2017 into branch

Location:
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90

    r8329 r9019  
    88   USE oce             ! ocean dynamics and tracers variables 
    99   USE dom_oce         ! ocean space and time domain 
     10   USE zdf_oce         ! ocean vertical physics     
     11   USE zdfgls   , ONLY : hmxl_n 
     12   ! 
    1013   USE in_out_manager  ! I/O units 
    1114   USE iom             ! I/0 library 
    12    USE wrk_nemo        ! working arrays 
    13 #if defined key_zdftke  
    14    USE zdf_oce, ONLY: en 
    15 #endif 
    16    USE zdf_oce, ONLY: avt, avm 
    17 #if defined key_zdfgls 
    18    USE zdf_oce, ONLY: en 
    19    USE zdfgls, ONLY: mxln 
    20 #endif 
    2115 
    2216   IMPLICIT NONE 
    2317   PRIVATE 
    2418 
    25    LOGICAL , PUBLIC ::   ln_dia25h     !:  25h mean output 
    2619   PUBLIC   dia_25h_init               ! routine called by nemogcm.F90 
    2720   PUBLIC   dia_25h                    ! routine called by diawri.F90 
    2821 
    29   !! * variables for calculating 25-hourly means 
    30    REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h 
    31    REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h  
    32    REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h 
    33    REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_25h , avm_25h 
    34 #if defined key_zdfgls || key_zdftke 
    35    REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_25h 
    36 #endif 
    37 #if defined key_zdfgls  
    38    REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   rmxln_25h 
    39 #endif 
    40    INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means 
    41  
    42  
     22   LOGICAL, PUBLIC ::   ln_dia25h      !:  25h mean output 
     23 
     24   ! variables for calculating 25-hourly means 
     25   INTEGER , SAVE ::   cnt_25h           ! Counter for 25 hour means 
     26   REAL(wp), SAVE ::   r1_25 = 0.04_wp   ! =1/25  
     27   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   tn_25h  , sn_25h 
     28   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   sshn_25h  
     29   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h 
     30   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avt_25h , avm_25h 
     31   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   en_25h  , rmxln_25h 
    4332 
    4433   !!---------------------------------------------------------------------- 
     
    5645      !!         
    5746      !! ** Method : Read namelist 
    58       !!   History 
    59       !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_25h 
    6047      !!--------------------------------------------------------------------------- 
    61       !! 
    6248      INTEGER ::   ios                 ! Local integer output status for namelist read 
    6349      INTEGER ::   ierror              ! Local integer for memory allocation 
     
    7965         WRITE(numout,*) 'dia_25h_init : Output 25 hour mean diagnostics' 
    8066         WRITE(numout,*) '~~~~~~~~~~~~' 
    81          WRITE(numout,*) 'Namelist nam_dia25h : set 25h outputs ' 
    82          WRITE(numout,*) 'Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h 
     67         WRITE(numout,*) '   Namelist nam_dia25h : set 25h outputs ' 
     68         WRITE(numout,*) '      Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h 
    8369      ENDIF 
    8470      IF( .NOT. ln_dia25h )   RETURN 
     
    8672      ! 1 - Allocate memory ! 
    8773      ! ------------------- ! 
    88       ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror ) 
     74      !                                ! ocean arrays 
     75      ALLOCATE( tn_25h (jpi,jpj,jpk), sn_25h (jpi,jpj,jpk), sshn_25h(jpi,jpj)  ,     & 
     76         &      un_25h (jpi,jpj,jpk), vn_25h (jpi,jpj,jpk), wn_25h(jpi,jpj,jpk),     & 
     77         &      avt_25h(jpi,jpj,jpk), avm_25h(jpi,jpj,jpk),                      STAT=ierror ) 
    8978      IF( ierror > 0 ) THEN 
    90          CALL ctl_stop( 'dia_25h: unable to allocate tn_25h' )   ;   RETURN 
    91       ENDIF 
    92       ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror ) 
    93       IF( ierror > 0 ) THEN 
    94          CALL ctl_stop( 'dia_25h: unable to allocate sn_25h' )   ;   RETURN 
    95       ENDIF 
    96       ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror ) 
    97       IF( ierror > 0 ) THEN 
    98          CALL ctl_stop( 'dia_25h: unable to allocate un_25h' )   ;   RETURN 
    99       ENDIF 
    100       ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror ) 
    101       IF( ierror > 0 ) THEN 
    102          CALL ctl_stop( 'dia_25h: unable to allocate vn_25h' )   ;   RETURN 
    103       ENDIF 
    104       ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror ) 
    105       IF( ierror > 0 ) THEN 
    106          CALL ctl_stop( 'dia_25h: unable to allocate wn_25h' )   ;   RETURN 
    107       ENDIF 
    108       ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror ) 
    109       IF( ierror > 0 ) THEN 
    110          CALL ctl_stop( 'dia_25h: unable to allocate avt_25h' )   ;   RETURN 
    111       ENDIF 
    112       ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror ) 
    113       IF( ierror > 0 ) THEN 
    114          CALL ctl_stop( 'dia_25h: unable to allocate avm_25h' )   ;   RETURN 
    115       ENDIF 
    116 # if defined key_zdfgls || defined key_zdftke 
    117       ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror ) 
    118       IF( ierror > 0 ) THEN 
    119          CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN 
    120       ENDIF 
    121 #endif 
    122 # if defined key_zdfgls  
    123       ALLOCATE( rmxln_25h(jpi,jpj,jpk), STAT=ierror ) 
    124       IF( ierror > 0 ) THEN 
    125          CALL ctl_stop( 'dia_25h: unable to allocate rmxln_25h' )   ;   RETURN 
    126       ENDIF 
    127 #endif 
    128       ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror ) 
    129       IF( ierror > 0 ) THEN 
    130          CALL ctl_stop( 'dia_25h: unable to allocate sshn_25h' )   ;   RETURN 
     79         CALL ctl_stop( 'dia_25h: unable to allocate ocean arrays' )   ;   RETURN 
     80      ENDIF 
     81      IF( ln_zdftke ) THEN             ! TKE physics 
     82         ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror ) 
     83         IF( ierror > 0 ) THEN 
     84            CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN 
     85         ENDIF 
     86      ENDIF 
     87      IF( ln_zdfgls ) THEN             ! GLS physics 
     88         ALLOCATE( en_25h(jpi,jpj,jpk), rmxln_25h(jpi,jpj,jpk), STAT=ierror ) 
     89         IF( ierror > 0 ) THEN 
     90            CALL ctl_stop( 'dia_25h: unable to allocate en_25h and rmxln_25h' )   ;   RETURN 
     91         ENDIF 
    13192      ENDIF 
    13293      ! ------------------------- ! 
     
    13495      ! ------------------------- ! 
    13596      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)  
    136       tn_25h(:,:,:) = tsb(:,:,:,jp_tem) 
    137       sn_25h(:,:,:) = tsb(:,:,:,jp_sal) 
    138       sshn_25h(:,:) = sshb(:,:) 
    139       un_25h(:,:,:) = ub(:,:,:) 
    140       vn_25h(:,:,:) = vb(:,:,:) 
    141       wn_25h(:,:,:) = wn(:,:,:) 
    142       avt_25h(:,:,:) = avt(:,:,:) 
    143       avm_25h(:,:,:) = avm(:,:,:) 
    144 # if defined key_zdfgls || defined key_zdftke 
     97      tn_25h  (:,:,:) = tsb (:,:,:,jp_tem) 
     98      sn_25h  (:,:,:) = tsb (:,:,:,jp_sal) 
     99      sshn_25h(:,:)   = sshb(:,:) 
     100      un_25h  (:,:,:) = ub  (:,:,:) 
     101      vn_25h  (:,:,:) = vb  (:,:,:) 
     102      wn_25h  (:,:,:) = wn  (:,:,:) 
     103      avt_25h (:,:,:) = avt (:,:,:) 
     104      avm_25h (:,:,:) = avm (:,:,:) 
     105      IF( ln_zdftke ) THEN 
    145106         en_25h(:,:,:) = en(:,:,:) 
    146 #endif 
    147 # if defined key_zdfgls 
    148          rmxln_25h(:,:,:) = mxln(:,:,:) 
    149 #endif 
    150 #if defined key_lim3 || defined key_lim2 
    151          CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice') 
     107      ENDIF 
     108      IF( ln_zdfgls ) THEN 
     109         en_25h   (:,:,:) = en    (:,:,:) 
     110         rmxln_25h(:,:,:) = hmxl_n(:,:,:) 
     111      ENDIF 
     112#if defined key_lim3 
     113      CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice') 
    152114#endif  
    153  
    154       ! -------------------------- ! 
    155       ! 3 - Return to dia_wri      ! 
    156       ! -------------------------- ! 
    157  
    158  
     115      ! 
    159116   END SUBROUTINE dia_25h_init 
    160117 
     
    164121      !!                 ***  ROUTINE dia_25h  *** 
    165122      !!          
    166       !! 
    167       !!-------------------------------------------------------------------- 
    168       !!                    
    169123      !! ** Purpose :   Write diagnostics with M2/S2 tide removed 
    170124      !! 
    171       !! ** Method  :    
    172       !!      25hr mean outputs for shelf seas 
     125      !! ** Method  :   25hr mean outputs for shelf seas 
     126      !!---------------------------------------------------------------------- 
     127      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    173128      !! 
    174       !! History : 
    175       !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam 
    176       !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam 
    177       !!   3.6  !  08-14  (E. O'Dea) adapted for VN3.6 
    178       !!---------------------------------------------------------------------- 
    179       !! * Modules used 
    180  
    181       IMPLICIT NONE 
    182  
    183       !! * Arguments 
    184       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    185  
    186  
    187       !! * Local declarations 
    188129      INTEGER ::   ji, jj, jk 
    189  
     130      INTEGER                          ::   iyear0, nimonth0,iday0            ! start year,imonth,day 
    190131      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout 
    191       REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi       ! temporary reals 
    192       INTEGER                          ::   i_steps                               ! no of timesteps per hour 
    193       REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace 
    194       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace 
    195       REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace 
    196       INTEGER                          ::   iyear0, nimonth0,iday0                ! start year,imonth,day 
    197  
     132      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! local scalars 
     133      INTEGER                          ::   i_steps                           ! no of timesteps per hour 
     134      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! workspace 
     135      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! workspace 
     136      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! workspace 
    198137      !!---------------------------------------------------------------------- 
    199138 
     
    207146      ENDIF 
    208147 
    209 #if defined key_lim3 || defined key_lim2 
    210       CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice') 
    211 #endif 
    212  
    213148      ! local variable for debugging 
    214149      ll_print = ll_print .AND. lwp 
    215150 
    216       ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours 
    217       ! every day 
    218       IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN 
     151      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours every day 
     152      IF( MOD( kt, i_steps ) == 0  .AND. kt /= nn_it000 ) THEN 
    219153 
    220154         IF (lwp) THEN 
     
    223157         ENDIF 
    224158 
    225          tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem) 
    226          sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal) 
    227          sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:) 
    228          un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:) 
    229          vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:) 
    230          wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:) 
    231          avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:) 
    232          avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:) 
    233 # if defined key_zdfgls || defined key_zdftke 
    234          en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:) 
    235 #endif 
    236 # if defined key_zdfgls 
    237          rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:) 
    238 #endif 
     159         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem) 
     160         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal) 
     161         sshn_25h(:,:)       = sshn_25h(:,:)   + sshn(:,:) 
     162         un_25h  (:,:,:)     = un_25h  (:,:,:) + un  (:,:,:) 
     163         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vn  (:,:,:) 
     164         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + wn  (:,:,:) 
     165         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:) 
     166         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:) 
     167         IF( ln_zdftke ) THEN 
     168            en_25h(:,:,:)    = en_25h  (:,:,:) + en(:,:,:) 
     169         ENDIF 
     170         IF( ln_zdfgls ) THEN 
     171            en_25h   (:,:,:) = en_25h   (:,:,:) + en    (:,:,:) 
     172            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxl_n(:,:,:) 
     173         ENDIF 
    239174         cnt_25h = cnt_25h + 1 
    240  
     175         ! 
    241176         IF (lwp) THEN 
    242177            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h 
    243178         ENDIF 
    244  
     179         ! 
    245180      ENDIF ! MOD( kt, i_steps ) == 0 
    246181 
    247          ! Write data for 25 hour mean output streams 
    248       IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 
    249  
    250             IF(lwp) THEN 
    251                WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt 
    252                WRITE(numout,*) '~~~~~~~~~~~~ ' 
    253             ENDIF 
    254  
    255             tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp 
    256             sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp 
    257             sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp 
    258             un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp 
    259             vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp 
    260             wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp 
    261             avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp 
    262             avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp 
    263 # if defined key_zdfgls || defined key_zdftke 
    264             en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp 
    265 #endif 
    266 # if defined key_zdfgls 
    267             rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp 
    268 #endif 
    269  
    270             IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output' 
    271             zmdi=1.e+20 !missing data indicator for masking 
    272             ! write tracers (instantaneous) 
    273             zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    274             CALL iom_put("temper25h", zw3d)   ! potential temperature 
    275             zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    276             CALL iom_put( "salin25h", zw3d  )   ! salinity 
    277             zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
    278             CALL iom_put( "ssh25h", zw2d )   ! sea surface  
    279  
    280  
    281             ! Write velocities (instantaneous) 
    282             zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:)) 
    283             CALL iom_put("vozocrtx25h", zw3d)    ! i-current 
    284             zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:)) 
    285             CALL iom_put("vomecrty25h", zw3d  )   ! j-current 
    286  
    287             zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    288             CALL iom_put("vomecrtz25h", zw3d )   ! k-current 
    289             zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    290             CALL iom_put("avt25h", zw3d )   ! diffusivity 
    291             zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    292             CALL iom_put("avm25h", zw3d)   ! viscosity 
    293 #if defined key_zdftke || defined key_zdfgls  
    294             zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     182      ! Write data for 25 hour mean output streams 
     183      IF( cnt_25h == 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt /= nn_it000 ) THEN 
     184         ! 
     185         IF(lwp) THEN 
     186            WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt 
     187            WRITE(numout,*) '~~~~~~~~~~~~ ' 
     188         ENDIF 
     189         ! 
     190         tn_25h  (:,:,:) = tn_25h  (:,:,:) * r1_25 
     191         sn_25h  (:,:,:) = sn_25h  (:,:,:) * r1_25 
     192         sshn_25h(:,:)   = sshn_25h(:,:)   * r1_25 
     193         un_25h  (:,:,:) = un_25h  (:,:,:) * r1_25 
     194         vn_25h  (:,:,:) = vn_25h  (:,:,:) * r1_25 
     195         wn_25h  (:,:,:) = wn_25h  (:,:,:) * r1_25 
     196         avt_25h (:,:,:) = avt_25h (:,:,:) * r1_25 
     197         avm_25h (:,:,:) = avm_25h (:,:,:) * r1_25 
     198         IF( ln_zdftke ) THEN 
     199            en_25h(:,:,:) = en_25h(:,:,:) * r1_25 
     200         ENDIF 
     201         IF( ln_zdfgls ) THEN 
     202            en_25h   (:,:,:) = en_25h   (:,:,:) * r1_25 
     203            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) * r1_25 
     204         ENDIF 
     205         ! 
     206         IF(lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output' 
     207         zmdi=1.e+20 !missing data indicator for masking 
     208         ! write tracers (instantaneous) 
     209         zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     210         CALL iom_put("temper25h", zw3d)   ! potential temperature 
     211         zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     212         CALL iom_put( "salin25h", zw3d  )   ! salinity 
     213         zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     214         CALL iom_put( "ssh25h", zw2d )   ! sea surface  
     215         ! Write velocities (instantaneous) 
     216         zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:)) 
     217         CALL iom_put("vozocrtx25h", zw3d)    ! i-current 
     218         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:)) 
     219         CALL iom_put("vomecrty25h", zw3d  )   ! j-current 
     220         zw3d(:,:,:) = wn_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     221         CALL iom_put("vomecrtz25h", zw3d )   ! k-current 
     222         ! Write vertical physics 
     223         zw3d(:,:,:) = avt_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     224         CALL iom_put("avt25h", zw3d )   ! diffusivity 
     225         zw3d(:,:,:) = avm_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     226         CALL iom_put("avm25h", zw3d)   ! viscosity 
     227         IF( ln_zdftke ) THEN 
     228            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    295229            CALL iom_put("tke25h", zw3d)   ! tke 
    296 #endif 
    297 #if defined key_zdfgls  
    298             zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     230         ENDIF 
     231         IF( ln_zdfgls ) THEN 
     232            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     233            CALL iom_put("tke25h", zw3d)   ! tke 
     234            zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    299235            CALL iom_put( "mxln25h",zw3d) 
    300 #endif 
    301  
    302             ! After the write reset the values to cnt=1 and sum values equal current value  
    303             tn_25h(:,:,:) = tsn(:,:,:,jp_tem) 
    304             sn_25h(:,:,:) = tsn(:,:,:,jp_sal) 
    305             sshn_25h(:,:) = sshn (:,:) 
    306             un_25h(:,:,:) = un(:,:,:) 
    307             vn_25h(:,:,:) = vn(:,:,:) 
    308             wn_25h(:,:,:) = wn(:,:,:) 
    309             avt_25h(:,:,:) = avt(:,:,:) 
    310             avm_25h(:,:,:) = avm(:,:,:) 
    311 # if defined key_zdfgls || defined key_zdftke 
     236         ENDIF 
     237         ! 
     238         ! After the write reset the values to cnt=1 and sum values equal current value  
     239         tn_25h  (:,:,:) = tsn (:,:,:,jp_tem) 
     240         sn_25h  (:,:,:) = tsn (:,:,:,jp_sal) 
     241         sshn_25h(:,:)   = sshn(:,:) 
     242         un_25h  (:,:,:) = un  (:,:,:) 
     243         vn_25h  (:,:,:) = vn  (:,:,:) 
     244         wn_25h  (:,:,:) = wn  (:,:,:) 
     245         avt_25h (:,:,:) = avt (:,:,:) 
     246         avm_25h (:,:,:) = avm (:,:,:) 
     247         IF( ln_zdftke ) THEN 
    312248            en_25h(:,:,:) = en(:,:,:) 
    313 #endif 
    314 # if defined key_zdfgls 
    315             rmxln_25h(:,:,:) = mxln(:,:,:) 
    316 #endif 
    317             cnt_25h = 1 
    318             IF (lwp)  WRITE(numout,*) 'dia_wri_tide :   & 
    319         &    After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h 
    320  
     249         ENDIF 
     250         IF( ln_zdfgls ) THEN 
     251            en_25h   (:,:,:) = en    (:,:,:) 
     252            rmxln_25h(:,:,:) = hmxl_n(:,:,:) 
     253         ENDIF 
     254         cnt_25h = 1 
     255         IF(lwp)  WRITE(numout,*) 'dia_wri_tide :   & 
     256            &    After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average', cnt_25h 
    321257      ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000 
    322  
    323  
     258      ! 
    324259   END SUBROUTINE dia_25h  
    325260 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r8083 r9019  
    3939       
    4040   !! * Substitutions 
    41 #  include "zdfddm_substitute.h90" 
    4241#  include "vectopt_loop_substitute.h90" 
    4342   !!---------------------------------------------------------------------- 
     
    214213         CALL wrk_alloc( jpi, jpj, zpe ) 
    215214         zpe(:,:) = 0._wp 
    216          IF( lk_zdfddm ) THEN 
     215         IF( ln_zdfddm ) THEN 
    217216            DO jk = 2, jpk 
    218217               DO jj = 1, jpj 
     
    221220                        zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   & 
    222221                           &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
     222!!gm  this can be reduced to :  (depw-dept) / e3w   (NB idem dans bn2 !) 
     223!                        zrw =   ( gdept_n(ji,jj,jk) - gdepw_n(ji,jj,jk) ) / e3w_n(ji,jj,jk) 
     224!!gm end 
    223225                        ! 
    224226                        zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw 
     
    226228                        ! 
    227229                        zpe(ji, jj) = zpe(ji, jj)            & 
    228                            &        -  grav * (    avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  & 
    229                            &                   - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) ) 
     230                           &        -  grav * (  avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  & 
     231                           &                   - avs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) ) 
    230232                     ENDIF 
    231233                  END DO 
     
    241243            END DO 
    242244         ENDIF 
    243          CALL lbc_lnk( zpe, 'T', 1._wp)          
    244          CALL iom_put( 'tnpeo', zpe ) 
    245          CALL wrk_dealloc( jpi, jpj, zpe ) 
     245!!gm useless lbc_lnk since the computation above is performed over 1:jpi & 1:jpj 
     246!!gm           CALL lbc_lnk( zpe, 'T', 1._wp)          
     247          CALL iom_put( 'tnpeo', zpe ) 
     248          CALL wrk_dealloc( jpi, jpj, zpe ) 
    246249      ENDIF 
    247250      ! 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90

    r8329 r9019  
    11MODULE diacfl 
    2    !!============================================================================== 
     2   !!====================================================================== 
    33   !!                       ***  MODULE  diacfl  *** 
    44   !! Output CFL diagnostics to ascii file 
    5    !!============================================================================== 
    6    !! History :  1.0  !  2010-03  (E. Blockley)  Original code 
    7    !!                 !  2014-06  (T Graham) Removed CPP key & Updated to vn3.6 
    8    !!  
     5   !!====================================================================== 
     6   !! History :  3.4  !  2010-03  (E. Blockley)  Original code 
     7   !!            3.6  !  2014-06  (T. Graham) Removed CPP key & Updated to vn3.6 
     8   !!            4.0  !  2017-09  (G. Madec)  style + comments 
    99   !!---------------------------------------------------------------------- 
    1010   !!   dia_cfl        : Compute and output Courant numbers at each timestep 
     
    1212   USE oce             ! ocean dynamics and active tracers 
    1313   USE dom_oce         ! ocean space and time domain 
     14   USE domvvl          !  
     15   ! 
    1416   USE lib_mpp         ! distribued memory computing 
    1517   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    1618   USE in_out_manager  ! I/O manager 
    17    USE domvvl      
    1819   USE timing          ! Performance output 
    1920 
     
    2122   PRIVATE 
    2223 
    23    REAL(wp) :: cu_max, cv_max, cw_max                      ! Run max U Courant number  
    24    INTEGER, DIMENSION(3) :: cu_loc, cv_loc, cw_loc         ! Run max locations 
    25    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcu_cfl           ! Courant number arrays 
    26    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcv_cfl           ! Courant number arrays 
    27    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcw_cfl           ! Courant number arrays 
     24   CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"    ! ascii filename 
     25   INTEGER           :: numcfl                            ! outfile unit 
     26   ! 
     27   INTEGER, DIMENSION(3) ::   nCu_loc, nCv_loc, nCw_loc   ! U, V, and W run max locations in the global domain 
     28   REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number  
    2829 
    29    INTEGER  :: numcfl                                       ! outfile unit 
    30    CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"      ! ascii filename 
     30!!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc ! 
     31!!gm          I don't understand why. 
     32   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
     33!!gm end 
    3134 
    3235   PUBLIC   dia_cfl       ! routine called by step.F90 
     
    4043   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4144   !!---------------------------------------------------------------------- 
    42  
    43  
    4445CONTAINS 
    45  
    4646 
    4747   SUBROUTINE dia_cfl ( kt ) 
     
    5252      !!               and output to ascii file 'cfl_diagnostics.ascii' 
    5353      !!---------------------------------------------------------------------- 
     54      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     55      ! 
     56      INTEGER :: ji, jj, jk   ! dummy loop indices 
     57      REAL(wp)::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars 
     58      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc   ! workspace 
     59!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
     60      !!---------------------------------------------------------------------- 
     61      ! 
     62      IF( nn_timing == 1 )   CALL timing_start('dia_cfl') 
     63      ! 
     64      !                       ! setup timestep multiplier to account for initial Eulerian timestep 
     65      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2dt = rdt 
     66      ELSE                                        ;    z2dt = rdt * 2._wp 
     67      ENDIF 
     68      ! 
     69      !                 
     70      DO jk = 1, jpk       ! calculate Courant numbers 
     71         DO jj = 1, jpj 
     72            DO ji = 1, fs_jpim1   ! vector opt. 
     73               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction 
     74               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction 
     75               zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk)   ! for k-direction 
     76            END DO 
     77         END DO          
     78      END DO 
     79      ! 
     80      !                    ! calculate maximum values and locations 
     81      IF( lk_mpp ) THEN 
     82         CALL mpp_maxloc( zCu_cfl, umask, zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) ) 
     83         CALL mpp_maxloc( zCv_cfl, vmask, zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) ) 
     84         CALL mpp_maxloc( zCw_cfl, wmask, zCw_max, iloc_w(1), iloc_w(2), iloc_w(3) ) 
     85      ELSE 
     86         iloc = MAXLOC( ABS( zcu_cfl(:,:,:) ) ) 
     87         iloc_u(1) = iloc(1) + nimpp - 1 
     88         iloc_u(2) = iloc(2) + njmpp - 1 
     89         iloc_u(3) = iloc(3) 
     90         zCu_max = zCu_cfl(iloc(1),iloc(2),iloc(3)) 
     91         ! 
     92         iloc = MAXLOC( ABS( zcv_cfl(:,:,:) ) ) 
     93         iloc_v(1) = iloc(1) + nimpp - 1 
     94         iloc_v(2) = iloc(2) + njmpp - 1 
     95         iloc_v(3) = iloc(3) 
     96         zCv_max = zCv_cfl(iloc(1),iloc(2),iloc(3)) 
     97         ! 
     98         iloc = MAXLOC( ABS( zcw_cfl(:,:,:) ) ) 
     99         iloc_w(1) = iloc(1) + nimpp - 1 
     100         iloc_w(2) = iloc(2) + njmpp - 1 
     101         iloc_w(3) = iloc(3) 
     102         zCw_max = zCw_cfl(iloc(1),iloc(2),iloc(3)) 
     103      ENDIF 
     104      ! 
     105      !                    ! write out to file 
     106      IF( lwp ) THEN 
     107         WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 
     108         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) 
     109         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3) 
     110      ENDIF 
     111      ! 
     112      !                    ! update maximum Courant numbers from whole run if applicable 
     113      IF( zCu_max > rCu_max ) THEN   ;   rCu_max = zCu_max   ;   nCu_loc(:) = iloc_u(:)   ;   ENDIF 
     114      IF( zCv_max > rCv_max ) THEN   ;   rCv_max = zCv_max   ;   nCv_loc(:) = iloc_v(:)   ;   ENDIF 
     115      IF( zCw_max > rCw_max ) THEN   ;   rCw_max = zCw_max   ;   nCw_loc(:) = iloc_w(:)   ;   ENDIF 
    54116 
    55       INTEGER, INTENT(in) ::  kt                            ! ocean time-step index 
     117      !                    ! at end of run output max Cu and Cv and close ascii file 
     118      IF( kt == nitend .AND. lwp ) THEN 
     119         ! to ascii file 
     120         WRITE(numcfl,*) '******************************************' 
     121         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3) 
     122         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCu_max 
     123         WRITE(numcfl,*) '******************************************' 
     124         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3) 
     125         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCv_max 
     126         WRITE(numcfl,*) '******************************************' 
     127         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3) 
     128         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCw_max 
     129         CLOSE( numcfl )  
     130         ! 
     131         ! to ocean output 
     132         WRITE(numout,*) 
     133         WRITE(numout,*) 'dia_cfl : Maximum Courant number information for the run ' 
     134         WRITE(numout,*) '~~~~~~~' 
     135         WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', z2dt/rCu_max 
     136         WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', z2dt/rCv_max 
     137         WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', z2dt/rCw_max 
     138      ENDIF 
     139      ! 
     140      IF( nn_timing == 1 )   CALL timing_stop('dia_cfl') 
     141      ! 
     142   END SUBROUTINE dia_cfl 
    56143 
    57       REAL(wp) :: zcu_max, zcv_max, zcw_max                 ! max Courant numbers per timestep 
    58       INTEGER, DIMENSION(3) :: zcu_loc, zcv_loc, zcw_loc    ! max Courant number locations 
    59  
    60       REAL(wp) :: dt                                        ! temporary scalars 
    61       INTEGER, DIMENSION(3) :: zlocu, zlocv, zlocw          ! temporary arrays  
    62       INTEGER  :: ji, jj, jk                                ! dummy loop indices 
    63  
    64        
    65       IF( nn_diacfl == 1) THEN 
    66          IF( nn_timing == 1 )   CALL timing_start('dia_cfl') 
    67          ! setup timestep multiplier to account for initial Eulerian timestep 
    68          IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    dt = rdt 
    69          ELSE                                        ;    dt = rdt * 2.0 
    70          ENDIF 
    71  
    72              ! calculate Courant numbers 
    73          DO jk = 1, jpk 
    74             DO jj = 1, jpj 
    75                DO ji = 1, fs_jpim1   ! vector opt. 
    76  
    77                   ! Courant number for x-direction (zonal current) 
    78                   zcu_cfl(ji,jj,jk) = ABS(un(ji,jj,jk))*dt/e1u(ji,jj) 
    79  
    80                   ! Courant number for y-direction (meridional current) 
    81                   zcv_cfl(ji,jj,jk) = ABS(vn(ji,jj,jk))*dt/e2v(ji,jj) 
    82  
    83                   ! Courant number for z-direction (vertical current) 
    84                   zcw_cfl(ji,jj,jk) = ABS(wn(ji,jj,jk))*dt/e3w_n(ji,jj,jk) 
    85                END DO 
    86             END DO          
    87          END DO 
    88  
    89          ! calculate maximum values and locations 
    90          IF( lk_mpp ) THEN 
    91             CALL mpp_maxloc(zcu_cfl,umask,zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3)) 
    92             CALL mpp_maxloc(zcv_cfl,vmask,zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3)) 
    93             CALL mpp_maxloc(zcw_cfl,tmask,zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3)) 
    94          ELSE 
    95             zlocu = MAXLOC( ABS( zcu_cfl(:,:,:) ) ) 
    96             zcu_loc(1) = zlocu(1) + nimpp - 1 
    97             zcu_loc(2) = zlocu(2) + njmpp - 1 
    98             zcu_loc(3) = zlocu(3) 
    99             zcu_max = zcu_cfl(zcu_loc(1),zcu_loc(2),zcu_loc(3)) 
    100  
    101             zlocv = MAXLOC( ABS( zcv_cfl(:,:,:) ) ) 
    102             zcv_loc(1) = zlocv(1) + nimpp - 1 
    103             zcv_loc(2) = zlocv(2) + njmpp - 1 
    104             zcv_loc(3) = zlocv(3) 
    105             zcv_max = zcv_cfl(zcv_loc(1),zcv_loc(2),zcv_loc(3)) 
    106  
    107             zlocw = MAXLOC( ABS( zcw_cfl(:,:,:) ) ) 
    108             zcw_loc(1) = zlocw(1) + nimpp - 1 
    109             zcw_loc(2) = zlocw(2) + njmpp - 1 
    110             zcw_loc(3) = zlocw(3) 
    111             zcw_max = zcw_cfl(zcw_loc(1),zcw_loc(2),zcw_loc(3)) 
    112          ENDIF 
    113        
    114          ! write out to file 
    115          IF( lwp ) THEN 
    116             WRITE(numcfl,FMT='(2x,i4,5x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3) 
    117             WRITE(numcfl,FMT='(11x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cv', zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3) 
    118             WRITE(numcfl,FMT='(11x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cw', zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3) 
    119          ENDIF 
    120  
    121          ! update maximum Courant numbers from whole run if applicable 
    122          IF( zcu_max > cu_max ) THEN 
    123             cu_max = zcu_max 
    124             cu_loc = zcu_loc 
    125          ENDIF 
    126          IF( zcv_max > cv_max ) THEN 
    127             cv_max = zcv_max 
    128             cv_loc = zcv_loc 
    129          ENDIF 
    130          IF( zcw_max > cw_max ) THEN 
    131             cw_max = zcw_max 
    132             cw_loc = zcw_loc 
    133          ENDIF 
    134  
    135          ! at end of run output max Cu and Cv and close ascii file 
    136          IF( kt == nitend .AND. lwp ) THEN 
    137             ! to ascii file 
    138             WRITE(numcfl,*) '******************************************' 
    139             WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', cu_max, cu_loc(1), cu_loc(2), cu_loc(3) 
    140             WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max) 
    141             WRITE(numcfl,*) '******************************************' 
    142             WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', cv_max, cv_loc(1), cv_loc(2), cv_loc(3) 
    143             WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max) 
    144             WRITE(numcfl,*) '******************************************' 
    145             WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', cw_max, cw_loc(1), cw_loc(2), cw_loc(3) 
    146             WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max) 
    147             CLOSE( numcfl )  
    148  
    149             ! to ocean output 
    150             WRITE(numout,*) 
    151             WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:' 
    152             WRITE(numout,*) '~~~~~~~~~~~~' 
    153             WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cu', cu_max, 'at (i, j, k) =   & 
    154                           &   (', cu_loc(1), cu_loc(2), cu_loc(3), ')' 
    155             WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max) 
    156             WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cv', cv_max, 'at (i, j, k) =   & 
    157                           &   (', cv_loc(1), cv_loc(2), cv_loc(3), ')' 
    158             WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max) 
    159             WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cw', cw_max, 'at (i, j, k) =   & 
    160                           &   (', cw_loc(1), cw_loc(2), cw_loc(3), ')' 
    161             WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max) 
    162  
    163          ENDIF 
    164  
    165          IF( nn_timing == 1 )   CALL timing_stop('dia_cfl') 
    166       ENDIF 
    167  
    168    END SUBROUTINE dia_cfl 
    169144 
    170145   SUBROUTINE dia_cfl_init 
     
    174149      !! ** Purpose :   create output file, initialise arrays 
    175150      !!---------------------------------------------------------------------- 
    176  
    177  
    178       IF( nn_diacfl == 1 ) THEN 
    179          IF( nn_timing == 1 )   CALL timing_start('dia_cfl_init') 
    180  
    181          cu_max=0.0 
    182          cv_max=0.0 
    183          cw_max=0.0 
    184  
    185          ALLOCATE( zcu_cfl(jpi, jpj, jpk), zcv_cfl(jpi, jpj, jpk), zcw_cfl(jpi, jpj, jpk) ) 
    186  
    187          zcu_cfl(:,:,:)=0.0 
    188          zcv_cfl(:,:,:)=0.0 
    189          zcw_cfl(:,:,:)=0.0 
    190  
    191          IF( lwp ) THEN 
    192             WRITE(numout,*) 
    193             WRITE(numout,*) 'dia_cfl     : Outputting CFL diagnostics to '//TRIM(clname) 
    194             WRITE(numout,*) '~~~~~~~~~~~~' 
    195             WRITE(numout,*) 
    196  
    197             ! create output ascii file 
    198             CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 ) 
    199             WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k' 
    200             WRITE(numcfl,*) '******************************************' 
    201          ENDIF 
    202  
    203          IF( nn_timing == 1 )   CALL timing_stop('dia_cfl_init') 
    204  
     151      ! 
     152      IF(lwp) THEN 
     153         WRITE(numout,*) 
     154         WRITE(numout,*) 'dia_cfl : Outputting CFL diagnostics to ',TRIM(clname), ' file' 
     155         WRITE(numout,*) '~~~~~~~' 
     156         WRITE(numout,*) 
     157         ! 
     158         ! create output ascii file 
     159         CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 ) 
     160         WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k' 
     161         WRITE(numcfl,*) '******************************************' 
    205162      ENDIF 
    206  
     163      ! 
     164      rCu_max = 0._wp 
     165      rCv_max = 0._wp 
     166      rCw_max = 0._wp 
     167      ! 
     168!!gm required to work 
     169      ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) ) 
     170!!gm end 
     171      !       
    207172   END SUBROUTINE dia_cfl_init 
    208173 
     174   !!====================================================================== 
    209175END MODULE diacfl 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    r8126 r9019  
    3232   USE dianam          ! build name of file 
    3333   USE lib_mpp         ! distributed memory computing library 
    34 #if defined key_lim2 
    35    USE ice_2 
    36 #endif 
    3734#if defined key_lim3 
    3835   USE ice 
     
    747744           END DO !end of loop on the level 
    748745 
    749 #if defined key_lim2 || defined key_lim3 
     746#if defined key_lim3 
    750747 
    751748           !ICE CASE     
     
    769766              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 
    770767 
    771 #if defined key_lim2    
    772               transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   &  
    773                                    (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &  
    774                                   *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &  
    775                                     hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) 
    776               transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &  
    777                                     (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) 
    778 #endif 
    779768#if defined key_lim3 
    780769              DO jl=1,jpl 
    781                  transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*     & 
    782                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * & 
    783                                   ( ht_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  & 
    784                                     ht_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) ) 
     770                 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*       & 
     771                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * & 
     772                                  ( h_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  & 
     773                                    h_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) ) 
    785774                                    
    786775                 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   & 
    787                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) 
     776                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) 
    788777              END DO 
    789778#endif 
     
    956945                 ENDIF ! end of test if point is in class  
    957946     
    958               ENDDO ! end of loop on the classes  
    959   
    960            ENDDO ! loop over jk  
    961   
    962 #if defined key_lim2 || defined key_lim3  
     947              END DO ! end of loop on the classes  
     948  
     949           END DO ! loop over jk  
     950  
     951#if defined key_lim3  
    963952  
    964953           !ICE CASE      
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90

    r6140 r9019  
    139139      DO jj = 1, jpj 
    140140         DO ji = 1, jpi 
    141             zztmp = bathy(ji,jj) 
     141            zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    142142            hth     (ji,jj) = zztmp 
    143143            zabs2   (ji,jj) = zztmp 
     
    150150         DO jj = 1, jpj 
    151151            DO ji = 1, jpi 
    152                zztmp = bathy(ji,jj) 
     152               zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    153153               zrho0_3(ji,jj) = zztmp 
    154154               zrho0_1(ji,jj) = zztmp 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

    r7753 r9019  
    397397      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport 
    398398      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 ) 
    399 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp ) 
     399902   IF( ios > 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp ) 
    400400      IF(lwm) WRITE ( numond, namptr ) 
    401401 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r8465 r9019  
    2525   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields 
    2626   !!---------------------------------------------------------------------- 
    27    USE oce             ! ocean dynamics and tracers  
    28    USE dom_oce         ! ocean space and time domain 
    29    USE dynadv, ONLY: ln_dynadv_vec 
    30    USE zdf_oce         ! ocean vertical physics 
    31    USE ldftra          ! lateral physics: eddy diffusivity coef. 
    32    USE ldfdyn          ! lateral physics: eddy viscosity   coef. 
    33    USE sbc_oce         ! Surface boundary condition: ocean fields 
    34    USE sbc_ice         ! Surface boundary condition: ice fields 
    35    USE icb_oce         ! Icebergs 
    36    USE icbdia          ! Iceberg budgets 
    37    USE sbcssr          ! restoring term toward SST/SSS climatology 
    38    USE phycst          ! physical constants 
    39    USE zdfmxl          ! mixed layer 
    40    USE dianam          ! build name of file (routine) 
    41    USE zdfddm          ! vertical  physics: double diffusion 
    42    USE diahth          ! thermocline diagnostics 
    43    USE wet_dry         ! wetting and drying 
    44    USE sbcwave         ! wave parameters 
     27   USE oce            ! ocean dynamics and tracers  
     28   USE dom_oce        ! ocean space and time domain 
     29   USE phycst         ! physical constants 
     30   USE dianam         ! build name of file (routine) 
     31   USE diahth         ! thermocline diagnostics 
     32   USE dynadv   , ONLY: ln_dynadv_vec 
     33   USE icb_oce        ! Icebergs 
     34   USE icbdia         ! Iceberg budgets 
     35   USE ldftra         ! lateral physics: eddy diffusivity coef. 
     36   USE ldfdyn         ! lateral physics: eddy viscosity   coef. 
     37   USE sbc_oce        ! Surface boundary condition: ocean fields 
     38   USE sbc_ice        ! Surface boundary condition: ice fields 
     39   USE sbcssr         ! restoring term toward SST/SSS climatology 
     40   USE sbcwave        ! wave parameters 
     41   USE wet_dry        ! wetting and drying 
     42   USE zdf_oce        ! ocean vertical physics 
     43   USE zdfdrg         ! ocean vertical physics: top/bottom friction 
     44   USE zdfmxl         ! mixed layer 
    4545   ! 
    46    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    47    USE in_out_manager  ! I/O manager 
    48    USE diatmb          ! Top,middle,bottom output 
    49    USE dia25h          ! 25h Mean output 
    50    USE iom 
    51    USE ioipsl 
    52  
    53 #if defined key_lim2 
    54    USE limwri_2  
    55 #elif defined key_lim3 
    56    USE limwri  
     46   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     47   USE in_out_manager ! I/O manager 
     48   USE diatmb         ! Top,middle,bottom output 
     49   USE dia25h         ! 25h Mean output 
     50   USE iom            !  
     51   USE ioipsl         !  
     52 
     53#if defined key_lim3 
     54   USE icewri  
    5755#endif 
    5856   USE lib_mpp         ! MPP library 
     
    6058   USE diurnal_bulk    ! diurnal warm layer 
    6159   USE cool_skin       ! Cool skin 
    62    USE wrk_nemo        ! working array 
    6360 
    6461   IMPLICIT NONE 
     
    8077 
    8178   !! * Substitutions 
    82 #  include "zdfddm_substitute.h90" 
    8379#  include "vectopt_loop_substitute.h90" 
    8480   !!---------------------------------------------------------------------- 
     
    120116      !! ** Method  :  use iom_put 
    121117      !!---------------------------------------------------------------------- 
    122       !! 
    123118      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    124119      !! 
    125       INTEGER                      ::   ji, jj, jk              ! dummy loop indices 
    126       INTEGER                      ::   jkbot                   ! 
    127       REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !  
    128       !! 
    129       REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace 
    130       REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace 
     120      INTEGER ::   ji, jj, jk       ! dummy loop indices 
     121      INTEGER ::   ikbot            ! local integer 
     122      REAL(wp)::   zztmp , zztmpx   ! local scalar 
     123      REAL(wp)::   zztmp2, zztmpy   !   -      - 
     124      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace 
     125      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace 
    131126      !!---------------------------------------------------------------------- 
    132127      !  
    133128      IF( nn_timing == 1 )   CALL timing_start('dia_wri') 
    134129      !  
    135       CALL wrk_alloc( jpi , jpj      , z2d ) 
    136       CALL wrk_alloc( jpi , jpj, jpk , z3d ) 
    137       ! 
    138130      ! Output the initial state and forcings 
    139131      IF( ninist == 1 ) THEN                        
     
    163155         DO jj = 1, jpj 
    164156            DO ji = 1, jpi 
    165                jkbot = mbkt(ji,jj) 
    166                z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem) 
     157               ikbot = mbkt(ji,jj) 
     158               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem) 
    167159            END DO 
    168160         END DO 
     
    175167         DO jj = 1, jpj 
    176168            DO ji = 1, jpi 
    177                jkbot = mbkt(ji,jj) 
    178                z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal) 
     169               ikbot = mbkt(ji,jj) 
     170               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal) 
    179171            END DO 
    180172         END DO 
     
    183175 
    184176      IF ( iom_use("taubot") ) THEN                ! bottom stress 
     177         zztmp = rau0 * 0.25 
    185178         z2d(:,:) = 0._wp 
    186179         DO jj = 2, jpjm1 
    187180            DO ji = fs_2, fs_jpim1   ! vector opt. 
    188                zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  & 
    189                       &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )       
    190                zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  & 
    191                       &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  )  
    192                z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)  
     181               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   & 
     182                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   & 
     183                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   & 
     184                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2 
     185               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
    193186               ! 
    194             ENDDO 
    195          ENDDO 
     187            END DO 
     188         END DO 
    196189         CALL lbc_lnk( z2d, 'T', 1. ) 
    197190         CALL iom_put( "taubot", z2d )            
    198191      ENDIF 
    199192          
    200       CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current 
    201       CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current 
     193      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current 
     194      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current 
    202195      IF ( iom_use("sbu") ) THEN 
    203196         DO jj = 1, jpj 
    204197            DO ji = 1, jpi 
    205                jkbot = mbku(ji,jj) 
    206                z2d(ji,jj) = un(ji,jj,jkbot) 
     198               ikbot = mbku(ji,jj) 
     199               z2d(ji,jj) = un(ji,jj,ikbot) 
    207200            END DO 
    208201         END DO 
     
    210203      ENDIF 
    211204       
    212       CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current 
    213       CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current 
     205      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current 
     206      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current 
    214207      IF ( iom_use("sbv") ) THEN 
    215208         DO jj = 1, jpj 
    216209            DO ji = 1, jpi 
    217                jkbot = mbkv(ji,jj) 
    218                z2d(ji,jj) = vn(ji,jj,jkbot) 
     210               ikbot = mbkv(ji,jj) 
     211               z2d(ji,jj) = vn(ji,jj,ikbot) 
    219212            END DO 
    220213         END DO 
     
    233226      ENDIF 
    234227 
    235       CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef. 
    236       CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef. 
    237       CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm) 
    238  
    239       IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt  (:,:,:) ) ) ) 
    240       IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) ) 
     228      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef. 
     229      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef. 
     230      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef. 
     231 
     232      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) ) 
     233      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 
    241234 
    242235      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
     
    251244         END DO 
    252245         CALL lbc_lnk( z2d, 'T', 1. ) 
    253          CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient 
     246         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
    254247         z2d(:,:) = SQRT( z2d(:,:) ) 
    255          CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient 
     248         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient 
    256249      ENDIF 
    257250          
    258       ! clem: heat and salt content 
     251      ! heat and salt contents 
    259252      IF( iom_use("heatc") ) THEN 
    260253         z2d(:,:)  = 0._wp  
     
    266259            END DO 
    267260         END DO 
    268          CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2) 
     261         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
    269262      ENDIF 
    270263 
     
    278271            END DO 
    279272         END DO 
    280          CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2) 
     273         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    281274      ENDIF 
    282275      ! 
    283276      IF ( iom_use("eken") ) THEN 
    284          rke(:,:,jpk) = 0._wp                               !      kinetic energy  
     277         z3d(:,:,jk) = 0._wp  
    285278         DO jk = 1, jpkm1 
    286279            DO jj = 2, jpjm1 
    287280               DO ji = fs_2, fs_jpim1   ! vector opt. 
    288                   zztmp   = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    289                   zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e1e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)    & 
    290                      &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e1e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) )  & 
    291                      &          *  zztmp  
    292                   ! 
    293                   zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1e2v(ji,jj-1) * e3v_n(ji,jj-1,jk)    & 
    294                      &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1e2v(ji,jj  ) * e3v_n(ji,jj  ,jk) )  & 
    295                      &          *  zztmp  
    296                   ! 
    297                   rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy ) 
    298                   ! 
    299                ENDDO 
    300             ENDDO 
    301          ENDDO 
    302          CALL lbc_lnk( rke, 'T', 1. ) 
    303          CALL iom_put( "eken", rke )            
     281                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     282                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   & 
     283                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   & 
     284                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   & 
     285                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   ) 
     286               END DO 
     287            END DO 
     288         END DO 
     289         CALL lbc_lnk( z3d, 'T', 1. ) 
     290         CALL iom_put( "eken", z3d )                 ! kinetic energy 
    304291      ENDIF 
    305292      ! 
     
    313300            z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 
    314301         END DO 
    315          CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction 
    316          CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum 
     302         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction 
     303         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum 
    317304      ENDIF 
    318305       
    319306      IF( iom_use("u_heattr") ) THEN 
    320          z2d(:,:) = 0.e0  
     307         z2d(:,:) = 0._wp  
    321308         DO jk = 1, jpkm1 
    322309            DO jj = 2, jpjm1 
     
    327314         END DO 
    328315         CALL lbc_lnk( z2d, 'U', -1. ) 
    329          CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction 
     316         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction 
    330317      ENDIF 
    331318 
     
    340327         END DO 
    341328         CALL lbc_lnk( z2d, 'U', -1. ) 
    342          CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction 
     329         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction 
    343330      ENDIF 
    344331 
     
    349336            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 
    350337         END DO 
    351          CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction 
     338         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction 
    352339      ENDIF 
    353340       
     
    362349         END DO 
    363350         CALL lbc_lnk( z2d, 'V', -1. ) 
    364          CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction 
     351         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction 
    365352      ENDIF 
    366353 
    367354      IF( iom_use("v_salttr") ) THEN 
    368          z2d(:,:) = 0.e0  
     355         z2d(:,:) = 0._wp  
    369356         DO jk = 1, jpkm1 
    370357            DO jj = 2, jpjm1 
     
    375362         END DO 
    376363         CALL lbc_lnk( z2d, 'V', -1. ) 
    377          CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction 
    378       ENDIF 
    379  
    380       ! Vertical integral of temperature 
     364         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction 
     365      ENDIF 
     366 
    381367      IF( iom_use("tosmint") ) THEN 
    382          z2d(:,:)=0._wp 
     368         z2d(:,:) = 0._wp 
    383369         DO jk = 1, jpkm1 
    384370            DO jj = 2, jpjm1 
    385371               DO ji = fs_2, fs_jpim1   ! vector opt. 
    386                   z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem) 
     372                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem) 
    387373               END DO 
    388374            END DO 
    389375         END DO 
    390376         CALL lbc_lnk( z2d, 'T', -1. ) 
    391          CALL iom_put( "tosmint", z2d )  
    392       ENDIF 
    393  
    394       ! Vertical integral of salinity 
     377         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature 
     378      ENDIF 
    395379      IF( iom_use("somint") ) THEN 
    396380         z2d(:,:)=0._wp 
     
    398382            DO jj = 2, jpjm1 
    399383               DO ji = fs_2, fs_jpim1   ! vector opt. 
    400                   z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
     384                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
    401385               END DO 
    402386            END DO 
    403387         END DO 
    404388         CALL lbc_lnk( z2d, 'T', -1. ) 
    405          CALL iom_put( "somint", z2d )  
    406       ENDIF 
    407  
    408       CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2) 
    409       ! 
    410       CALL wrk_dealloc( jpi , jpj      , z2d ) 
    411       CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 
    412       ! 
    413       ! If we want tmb values  
    414  
    415       IF (ln_diatmb) THEN 
    416          CALL dia_tmb  
    417       ENDIF  
    418       IF (ln_dia25h) THEN 
    419          CALL dia_25h( kt ) 
    420       ENDIF  
     389         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity 
     390      ENDIF 
     391 
     392      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2) 
     393      ! 
     394 
     395      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values  
     396           
     397      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging 
    421398 
    422399      IF( nn_timing == 1 )   CALL timing_stop('dia_wri') 
     
    452429      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars 
    453430      ! 
    454       REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace 
    455       REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace 
     431      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace 
     432      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
    456433      !!---------------------------------------------------------------------- 
    457434      !  
    458435      IF( nn_timing == 1 )   CALL timing_start('dia_wri') 
    459436      ! 
    460                              CALL wrk_alloc( jpi,jpj      , zw2d ) 
    461       IF( .NOT.ln_linssh )   CALL wrk_alloc( jpi,jpj,jpk  , zw3d ) 
    462       ! 
    463       ! Output the initial state and forcings 
    464       IF( ninist == 1 ) THEN                        
     437      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==! 
    465438         CALL dia_wri_state( 'output.init', kt ) 
    466439         ninist = 0 
     
    470443      ! ----------------- 
    471444 
    472       ! local variable for debugging 
    473       ll_print = .FALSE. 
     445      ll_print = .FALSE.                  ! local variable for debugging 
    474446      ll_print = ll_print .AND. lwp 
    475447 
     
    707679#endif 
    708680 
    709          IF( ln_cpl .AND. nn_ice == 2 ) THEN 
    710             CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice 
    711                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    712             CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice 
    713                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    714          ENDIF 
    715  
    716681         CALL histend( nid_T, snc4chunks=snc4set ) 
    717682 
     
    747712         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt 
    748713            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
    749          CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu 
     714         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm 
    750715            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
    751716 
    752          IF( lk_zdfddm ) THEN 
     717         IF( ln_zdfddm ) THEN 
    753718            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs 
    754719               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
     
    861826#endif 
    862827 
    863       IF( ln_cpl .AND. nn_ice == 2 ) THEN 
    864          CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature 
    865          CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo 
    866       ENDIF 
    867  
    868828      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current 
    869829      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress 
     
    874834      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current 
    875835      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef. 
    876       CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
    877       IF( lk_zdfddm ) THEN 
    878          CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef. 
     836      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
     837      IF( ln_zdfddm ) THEN 
     838         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef. 
    879839      ENDIF 
    880840 
    881841      IF( ln_wave .AND. ln_sdw ) THEN 
    882          CALL histwrite( nid_U, "sdzocrtx", it, usd           , ndim_U , ndex_U )    ! i-StokesDrift-current 
    883          CALL histwrite( nid_V, "sdmecrty", it, vsd           , ndim_V , ndex_V )    ! j-StokesDrift-current 
    884          CALL histwrite( nid_W, "sdvecrtz", it, wsd           , ndim_T , ndex_T )    ! StokesDrift vert. current 
     842         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current 
     843         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current 
     844         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current 
    885845      ENDIF 
    886846 
     
    893853         CALL histclo( nid_W ) 
    894854      ENDIF 
    895       ! 
    896                              CALL wrk_dealloc( jpi , jpj        , zw2d ) 
    897       IF( .NOT.ln_linssh )   CALL wrk_dealloc( jpi , jpj , jpk  , zw3d ) 
    898855      ! 
    899856      IF( nn_timing == 1 )   CALL timing_stop('dia_wri') 
     
    1009966      ENDIF 
    1010967 
    1011 #if defined key_lim2 
    1012       CALL lim_wri_state_2( kt, id_i, nh_i ) 
    1013 #elif defined key_lim3 
    1014       CALL lim_wri_state( kt, id_i, nh_i ) 
     968#if defined key_lim3 
     969      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid 
     970         CALL ice_wri_state( kt, id_i, nh_i ) 
     971      ENDIF 
    1015972#else 
    1016973      CALL histend( id_i, snc4chunks=snc4set ) 
Note: See TracChangeset for help on using the changeset viewer.