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.
diaopfoam.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90 @ 10390

Last change on this file since 10390 was 10390, checked in by dford, 6 years ago

Add 25-hour mean and top-middle-bottom diagnostics for all variables available through FABM, and set up indices for specific variables for easier access by OBS/ASM. See internal Met Office NEMO ticket 760.
File size: 7.9 KB
Line 
1MODULE diaopfoam
2   !!======================================================================
3   !!                       ***  MODULE  diaopfoam  ***
4   !! Output stream for operational use
5   !!======================================================================
6   !! History :  3.6  !  2016  (P Sykes)  Original code
7   !!----------------------------------------------------------------------
8   USE oce             ! ocean dynamics and tracers variables
9   USE dom_oce         ! ocean space and time domain
10   USE diainsitutem, ONLY: rinsitu_t, theta2t
11   USE in_out_manager  ! I/O units
12   USE iom             ! I/0 library
13   USE wrk_nemo        ! working arrays
14   USE diatmb
15   USE diurnal_bulk
16   USE cool_skin
17
18
19   IMPLICIT NONE
20   PRIVATE
21
22   LOGICAL , PUBLIC ::   ln_diaopfoam     !: Diaopfoam output
23   PUBLIC   dia_diaopfoam_init            ! routine called by nemogcm.F90
24   PUBLIC   dia_diaopfoam                 ! routine called by diawri.F90
25   PUBLIC   calc_max_cur                  ! routine called by diaopfoam.F90
26
27   !! * Substitutions
28#  include "domzgr_substitute.h90"
29
30   !!----------------------------------------------------------------------
31   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
32   !! $Id$
33   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE dia_diaopfoam_init
38      !!---------------------------------------------------------------------------
39      !!                  ***  ROUTINE dia_wri_diaop_init  ***
40      !!     
41      !! ** Purpose: Initialization of diaopfoam namelist
42      !!       
43      !! ** Method : Read namelist
44      !!   History
45      !!   3.4  !  03-14  (P. Sykes) Routine to initialize dia_wri_diaop
46      !!---------------------------------------------------------------------------
47      !!
48      INTEGER            ::   ios      ! Local integer output status for namelist read
49      INTEGER            ::   ierror   ! local integer
50      !!
51      NAMELIST/nam_diadiaop/ ln_diaopfoam
52      !!
53      !!----------------------------------------------------------------------
54      !
55      ln_diaopfoam = .false.         ! default value for diaopfoam stream
56      REWIND ( numnam_ref )              ! Read Namelist nam_diadiaop in reference namelist : 3D hourly diagnostics
57      READ   ( numnam_ref, nam_diadiaop, IOSTAT=ios, ERR= 901 )
58901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in reference namelist', lwp )
59
60      REWIND( numnam_cfg )              ! Namelist nam_diadiaop in configuration namelist  3D hourly diagnostics
61      READ  ( numnam_cfg, nam_diadiaop, IOSTAT = ios, ERR = 902 )
62902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in configuration namelist', lwp )
63      IF(lwm) WRITE ( numond, nam_diadiaop )
64      !
65      IF(lwp) THEN                   ! Control print
66         WRITE(numout,*)
67         WRITE(numout,*) 'dia_diaopfoam_init : Output Diaopfoam Diagnostics'
68         WRITE(numout,*) '~~~~~~~~~~~~'
69         WRITE(numout,*) '   Namelist nam_diadiaop : set diaopfoam outputs '
70         WRITE(numout,*) '      Switch for diaopfoam diagnostics (T) or not (F)  ln_diaopfoam  = ', ln_diaopfoam
71      ENDIF
72    END SUBROUTINE dia_diaopfoam_init
73
74    SUBROUTINE dia_diaopfoam
75      !!----------------------------------------------------------------------
76      !!                 ***  ROUTINE dia_diaopfoam  ***
77      !! ** Purpose :   Write 3D hourly diagnostics for operational use
78      !!
79      !!
80      !! History :
81      !!   3.6  !  11-16  (P. Sykes)
82      !!         
83      !!--------------------------------------------------------------------
84      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
85      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
86      REAL(wp) :: zmdi
87      REAL(wp), POINTER, DIMENSION(:,:)   :: zwu
88      REAL(wp), POINTER, DIMENSION(:,:)   :: zwv
89      REAL(wp), POINTER, DIMENSION(:,:)   :: zwz
90
91      CALL wrk_alloc( jpi , jpj      , zwu )
92      CALL wrk_alloc( jpi , jpj      , zwv )
93      CALL wrk_alloc( jpi , jpj      , zwz )
94
95      zmdi=1.e+20                               !  missing data indicator for masking
96      ! Diaopfoam stream if needed
97      IF (ln_diaopfoam) THEN
98         CALL theta2t
99         CALL iom_put( "insitut_op" , rinsitu_t(:,:,:)                      )    ! insitu temperature
100         CALL iom_put( "toce_op"   , tsn(:,:,:,jp_tem)                     )    ! temperature
101         CALL iom_put( "soce_op"   , tsn(:,:,:,jp_sal)                     )    ! salinity
102         IF (ln_diurnal) THEN
103            CALL iom_put( "sst_wl_op"   , x_dsst                    )    ! warm layer
104            CALL iom_put( "sst_cs_op"   , x_csdsst                  )    ! cool skin
105         ENDIF
106         zw2d(:,:)=sshn(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
107         CALL iom_put( "ssh_op"  , zw2d(:,:)                          ) ! sea surface height
108         CALL iom_put( "uoce_op"   , un                                    )    ! i-current     
109         CALL iom_put( "voce_op"   , vn                                    )    ! j-current
110         !CALL iom_put( "woce_op"   , wn                                    )    ! k-current
111         CALL calc_max_cur(zwu,zwv,zwz,zmdi)
112         CALL iom_put( "maxu" , zwu                                     ) ! max u current
113         CALL iom_put( "maxv" , zwv                                     ) ! max v current
114         CALL iom_put( "maxz" , zwz                                     ) ! max current depth
115      ENDIF
116   END SUBROUTINE dia_diaopfoam
117
118   SUBROUTINE calc_max_cur(zmax_u, zmax_v, zmax_z, inmdi)
119      !!---------------------------------------------------------------------
120      !!                    ***  ROUTINE calc_max_cur  ***
121      !!
122      !! ** Purpose :   To locate within the water column the magnitude and
123      !!                vertical location of the strongest horizontal
124      !!                current. The vertical component is ignored since it
125      !!                is an order of magnitude smaller than the horizontal
126      !!                flow, in general.
127      !!
128      !! ** Method :    A. Map U,V to T-grid.
129      !!                B. Calculate the magnitude of the current for every
130      !!                   grid cell.
131      !!                C. Locate the vertical index using FORTRAN's builtin
132      !!                   MAXLOC function.
133      !!                D. Copy the U,V,Z components of the relevant
134      !!                   indices.
135      !!
136      !! ** Returns :   Value of u, v component and depth of maximum
137      !!                horizontal current on T-grid.
138      !!
139      !!---------------------------------------------------------------------
140      IMPLICIT NONE
141      REAL(wp), DIMENSION(jpi, jpj), INTENT(  OUT) :: zmax_u, zmax_v, zmax_z
142      REAL(wp), DIMENSION(jpi, jpj, jpk)           :: zmax_u_t, zmax_v_t
143      REAL(wp), DIMENSION(jpi, jpj, jpk)           :: zcmag
144      REAL(wp), DIMENSION(jpi, jpj)                :: zmaxk
145      REAL(wp),                      INTENT(IN   ) :: inmdi
146      INTEGER :: ji, jj, jk
147
148      ! Initialise output arrays
149      zmax_u(:, :) = inmdi
150      zmax_v(:, :) = inmdi
151      zmax_z(:, :) = inmdi
152      ! Map to T-grid
153      zmax_u_t(:, :, :) = 0._wp
154      zmax_v_t(:, :, :) = 0._wp
155      DO jk = 1,jpk
156         DO jj = 2,jpj
157            DO ji = 2,jpi
158               zmax_u_t(ji, jj, jk) = 0.5 * (un(ji, jj, jk) + un(ji-1, jj, jk))
159               zmax_v_t(ji, jj, jk) = 0.5 * (vn(ji, jj, jk) + vn(ji, jj-1, jk))
160            END DO
161         END DO
162      END DO
163      ! Calculate absolute velocity
164      zcmag = sqrt((zmax_u_t)**2 + (zmax_v_t)**2)
165      ! Find max. current
166      zmaxk = maxloc(zcmag, dim=3)
167      ! Output values
168      DO jj = 1,jpj
169         DO ji = 1,jpi
170             zmax_u(ji, jj) = zmax_u_t(ji, jj, INT(zmaxk(ji, jj)))
171             zmax_v(ji, jj) = zmax_v_t(ji, jj, INT(zmaxk(ji, jj)))
172             zmax_z(ji, jj) = fsdept(ji, jj, INT(zmaxk(ji, jj)))
173         END DO
174      END DO     
175   END SUBROUTINE calc_max_cur
176
177END MODULE diaopfoam
Note: See TracBrowser for help on using the repository browser.