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_bgc_updates/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_bgc_updates/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90 @ 10241

Last change on this file since 10241 was 10241, checked in by dford, 5 years ago

Add visibility diagnostic (can't currently be done via fabm.yaml without changing the FABM code).
File size: 8.1 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!#if defined key_spm
112!         cltra = TRIM(ctrc3d(5))//"_op"
113!         zw3d(:,:,:) = trc3d(:,:,:,5)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) ! Visibility
114!         CALL iom_put( cltra, zw3d  )
115!#endif
116         CALL calc_max_cur(zwu,zwv,zwz,zmdi)
117         CALL iom_put( "maxu" , zwu                                     ) ! max u current
118         CALL iom_put( "maxv" , zwv                                     ) ! max v current
119         CALL iom_put( "maxz" , zwz                                     ) ! max current depth
120      ENDIF
121   END SUBROUTINE dia_diaopfoam
122
123   SUBROUTINE calc_max_cur(zmax_u, zmax_v, zmax_z, inmdi)
124      !!---------------------------------------------------------------------
125      !!                    ***  ROUTINE calc_max_cur  ***
126      !!
127      !! ** Purpose :   To locate within the water column the magnitude and
128      !!                vertical location of the strongest horizontal
129      !!                current. The vertical component is ignored since it
130      !!                is an order of magnitude smaller than the horizontal
131      !!                flow, in general.
132      !!
133      !! ** Method :    A. Map U,V to T-grid.
134      !!                B. Calculate the magnitude of the current for every
135      !!                   grid cell.
136      !!                C. Locate the vertical index using FORTRAN's builtin
137      !!                   MAXLOC function.
138      !!                D. Copy the U,V,Z components of the relevant
139      !!                   indices.
140      !!
141      !! ** Returns :   Value of u, v component and depth of maximum
142      !!                horizontal current on T-grid.
143      !!
144      !!---------------------------------------------------------------------
145      IMPLICIT NONE
146      REAL(wp), DIMENSION(jpi, jpj), INTENT(  OUT) :: zmax_u, zmax_v, zmax_z
147      REAL(wp), DIMENSION(jpi, jpj, jpk)           :: zmax_u_t, zmax_v_t
148      REAL(wp), DIMENSION(jpi, jpj, jpk)           :: zcmag
149      REAL(wp), DIMENSION(jpi, jpj)                :: zmaxk
150      REAL(wp),                      INTENT(IN   ) :: inmdi
151      INTEGER :: ji, jj, jk
152
153      ! Initialise output arrays
154      zmax_u(:, :) = inmdi
155      zmax_v(:, :) = inmdi
156      zmax_z(:, :) = inmdi
157      ! Map to T-grid
158      zmax_u_t(:, :, :) = 0._wp
159      zmax_v_t(:, :, :) = 0._wp
160      DO jk = 1,jpk
161         DO jj = 2,jpj
162            DO ji = 2,jpi
163               zmax_u_t(ji, jj, jk) = 0.5 * (un(ji, jj, jk) + un(ji-1, jj, jk))
164               zmax_v_t(ji, jj, jk) = 0.5 * (vn(ji, jj, jk) + vn(ji, jj-1, jk))
165            END DO
166         END DO
167      END DO
168      ! Calculate absolute velocity
169      zcmag = sqrt((zmax_u_t)**2 + (zmax_v_t)**2)
170      ! Find max. current
171      zmaxk = maxloc(zcmag, dim=3)
172      ! Output values
173      DO jj = 1,jpj
174         DO ji = 1,jpi
175             zmax_u(ji, jj) = zmax_u_t(ji, jj, INT(zmaxk(ji, jj)))
176             zmax_v(ji, jj) = zmax_v_t(ji, jj, INT(zmaxk(ji, jj)))
177             zmax_z(ji, jj) = fsdept(ji, jj, INT(zmaxk(ji, jj)))
178         END DO
179      END DO     
180   END SUBROUTINE calc_max_cur
181
182END MODULE diaopfoam
Note: See TracBrowser for help on using the repository browser.