MODULE diaopfoam !!====================================================================== !! *** MODULE diaopfoam *** !! Output stream for operational use !!====================================================================== !! History : 3.6 ! 2016 (P Sykes) Original code !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain USE diainsitutem, ONLY: rinsitu_t, theta2t USE in_out_manager ! I/O units USE iom ! I/0 library USE wrk_nemo ! working arrays USE diatmb USE diurnal_bulk USE cool_skin USE ioipsl IMPLICIT NONE PRIVATE LOGICAL , PUBLIC :: ln_diaopfoam !: Diaopfoam output LOGICAL , PUBLIC :: ln_diaopfoam_Tzero !: Diaopfoam first time step output PUBLIC dia_diaopfoam_init ! routine called by nemogcm.F90 PUBLIC dia_diaopfoam ! routine called by diawri.F90 PUBLIC calc_max_cur ! routine called by diaopfoam.F90 !! * Substitutions # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.6 , NEMO Consortium (2014) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dia_diaopfoam_init !!--------------------------------------------------------------------------- !! *** ROUTINE dia_wri_diaop_init *** !! !! ** Purpose: Initialization of diaopfoam namelist !! !! ** Method : Read namelist !! History !! 3.4 ! 03-14 (P. Sykes) Routine to initialize dia_wri_diaop !!--------------------------------------------------------------------------- !! INTEGER :: ios ! Local integer output status for namelist read INTEGER :: ierror ! local integer !! NAMELIST/nam_diadiaop/ ln_diaopfoam,ln_diaopfoam_Tzero !! !!---------------------------------------------------------------------- ! ln_diaopfoam = .false. ! default value for diaopfoam stream ln_diaopfoam_Tzero = .false. ! default value for diaopfoam Tzero stream REWIND ( numnam_ref ) ! Read Namelist nam_diadiaop in reference namelist : 3D hourly diagnostics READ ( numnam_ref, nam_diadiaop, IOSTAT=ios, ERR= 901 ) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in reference namelist', lwp ) REWIND( numnam_cfg ) ! Namelist nam_diadiaop in configuration namelist 3D hourly diagnostics READ ( numnam_cfg, nam_diadiaop, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in configuration namelist', lwp ) IF(lwm) WRITE ( numond, nam_diadiaop ) ! IF(lwp) THEN ! Control print WRITE(numout,*) WRITE(numout,*) 'dia_diaopfoam_init : Output Diaopfoam Diagnostics' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) ' Namelist nam_diadiaop : set diaopfoam outputs ' WRITE(numout,*) ' Switch for diaopfoam diagnostics (T) or not (F) ln_diaopfoam = ', ln_diaopfoam WRITE(numout,*) ' Switch for diaopfoam first timestep diagnostics (T) or not (F) ln_diaopfoam_Tzero = ', ln_diaopfoam_Tzero ENDIF END SUBROUTINE dia_diaopfoam_init SUBROUTINE dia_diaopfoam( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dia_diaopfoam *** !! ** Purpose : Write 3D hourly diagnostics for operational use !! !! !! History : !! 3.6 ! 11-16 (P. Sykes) !! !!-------------------------------------------------------------------- IMPLICIT NONE INTEGER, INTENT( in ) :: kt ! ocean time-step index REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace REAL(wp) :: zmdi REAL(wp), POINTER, DIMENSION(:,:) :: zwu REAL(wp), POINTER, DIMENSION(:,:) :: zwv REAL(wp), POINTER, DIMENSION(:,:) :: zwz CALL wrk_alloc( jpi , jpj , zwu ) CALL wrk_alloc( jpi , jpj , zwv ) CALL wrk_alloc( jpi , jpj , zwz ) zmdi=1.e+20 ! missing data indicator for masking ! Diaopfoam stream if needed IF (ln_diaopfoam) THEN IF ( kt .eq. nn_it000 .AND. ln_diaopfoam_Tzero ) THEN IF(lwp) WRITE(numout,*) 'diaopfoam: writing T0 at kt = ', kt CALL dia_diaopfoam_zero( 'Tzero', kt ) ENDIF CALL theta2t CALL iom_put( "insitut_op" , rinsitu_t(:,:,:) ) ! insitu temperature CALL iom_put( "toce_op" , tsn(:,:,:,jp_tem) ) ! temperature CALL iom_put( "soce_op" , tsn(:,:,:,jp_sal) ) ! salinity IF (ln_diurnal) THEN CALL iom_put( "sst_wl_op" , x_dsst ) ! warm layer CALL iom_put( "sst_cs_op" , x_csdsst ) ! cool skin ENDIF zw2d(:,:)=sshn(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) CALL iom_put( "ssh_op" , zw2d(:,:) ) ! sea surface height CALL iom_put( "uoce_op" , un ) ! i-current CALL iom_put( "voce_op" , vn ) ! j-current !CALL iom_put( "woce_op" , wn ) ! k-current CALL calc_max_cur(zwu,zwv,zwz,zmdi) CALL iom_put( "maxu" , zwu ) ! max u current CALL iom_put( "maxv" , zwv ) ! max v current CALL iom_put( "maxz" , zwz ) ! max current depth ENDIF END SUBROUTINE dia_diaopfoam SUBROUTINE calc_max_cur(zmax_u, zmax_v, zmax_z, inmdi) !!--------------------------------------------------------------------- !! *** ROUTINE calc_max_cur *** !! !! ** Purpose : To locate within the water column the magnitude and !! vertical location of the strongest horizontal !! current. The vertical component is ignored since it !! is an order of magnitude smaller than the horizontal !! flow, in general. !! !! ** Method : A. Map U,V to T-grid. !! B. Calculate the magnitude of the current for every !! grid cell. !! C. Locate the vertical index using FORTRAN's builtin !! MAXLOC function. !! D. Copy the U,V,Z components of the relevant !! indices. !! !! ** Returns : Value of u, v component and depth of maximum !! horizontal current on T-grid. !! !!--------------------------------------------------------------------- IMPLICIT NONE REAL(wp), DIMENSION(jpi, jpj), INTENT( OUT) :: zmax_u, zmax_v, zmax_z REAL(wp), DIMENSION(jpi, jpj, jpk) :: zmax_u_t, zmax_v_t REAL(wp), DIMENSION(jpi, jpj, jpk) :: zcmag REAL(wp), DIMENSION(jpi, jpj) :: zmaxk REAL(wp), INTENT(IN ) :: inmdi INTEGER :: ji, jj, jk ! Initialise output arrays zmax_u(:, :) = inmdi zmax_v(:, :) = inmdi zmax_z(:, :) = inmdi ! Map to T-grid zmax_u_t(:, :, :) = 0._wp zmax_v_t(:, :, :) = 0._wp DO jk = 1,jpk DO jj = 2,jpj DO ji = 2,jpi zmax_u_t(ji, jj, jk) = 0.5 * (un(ji, jj, jk) + un(ji-1, jj, jk)) zmax_v_t(ji, jj, jk) = 0.5 * (vn(ji, jj, jk) + vn(ji, jj-1, jk)) END DO END DO END DO ! Calculate absolute velocity zcmag = sqrt((zmax_u_t)**2 + (zmax_v_t)**2) ! Find max. current zmaxk = maxloc(zcmag, dim=3) ! Output values DO jj = 1,jpj DO ji = 1,jpi zmax_u(ji, jj) = zmax_u_t(ji, jj, INT(zmaxk(ji, jj))) zmax_v(ji, jj) = zmax_v_t(ji, jj, INT(zmaxk(ji, jj))) zmax_z(ji, jj) = fsdept(ji, jj, INT(zmaxk(ji, jj))) END DO END DO END SUBROUTINE calc_max_cur SUBROUTINE dia_diaopfoam_zero( cdfile_name, kt ) !!--------------------------------------------------------------------- !! *** ROUTINE dia_diaopfoam_zero *** !! !! ** Purpose : create a NetCDF file named cdfile_name which contains !! the instantaneous ocean state at the first time step. !! !! ** Method : NetCDF files using ioipsl !!---------------------------------------------------------------------- CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created INTEGER , INTENT( in ) :: kt ! ocean time-step index !! CHARACTER (len=32) :: clhstnam CHARACTER (len=40) :: clop INTEGER :: iimi, iima, ipk, ijmi, ijma ! local integers INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file INTEGER :: id_i , nz_i, nh_i INTEGER, DIMENSION(1) :: idex ! local workspace INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V INTEGER :: ierr INTEGER :: jkbot, jj, ji REAL(wp) :: zsto, zout, zmax, zjulian, zdt REAL(wp) :: zmdi REAL(wp), POINTER, DIMENSION(:,:) :: zwu REAL(wp), POINTER, DIMENSION(:,:) :: zwv REAL(wp), POINTER, DIMENSION(:,:) :: zwz REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace REAL(wp), DIMENSION(jpi,jpj) :: z2d !!---------------------------------------------------------------------- ! ----------------------------------------------------------------- ! 0. Allocations ! ----------------------------------------------------------------- ierr = 0 ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) , & & ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) , & & ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr ) IF( lk_mpp ) CALL mpp_sum( ierr ) IF( ierr /= 0 ) THEN CALL ctl_stop('dia_diaopfoam_zero: failed to allocate arrays') RETURN ENDIF CALL wrk_alloc( jpi , jpj , zwu ) CALL wrk_alloc( jpi , jpj , zwv ) CALL wrk_alloc( jpi , jpj , zwz ) zmdi=1.e+20 ! missing data indicator for masking ! ----------------------------------------------------------------- ! 1. Define NETCDF files and fields at beginning of first time step ! ----------------------------------------------------------------- ! Define indices of the horizontal output zoom and vertical limit storage iimi = 1 ; iima = jpi ijmi = 1 ; ijma = jpj ipk = jpk ! Define frequency of output and means zdt = rdt zsto = rdt clop = "inst(x)" ! no use of the mask value (require less cpu time) zout = rdt zmax = ( nitend - nit000 + 1 ) * zdt ! Compute julian date from starting date of the run CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) ! time axis zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment ! Define the T grid FILE ( nid_T ) clhstnam = TRIM(cdfile_name)//".grid_T" IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & & nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept & "m", ipk, gdept_1d, nz_T, "down" ) ! ! Index of ocean points CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface ! Define the U grid FILE ( nid_U ) clhstnam = TRIM(cdfile_name)//".grid_U" IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & & nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdepu & "m", ipk, gdept_1d, nz_U, "down" ) ! ! Index of ocean points CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U ) ! volume CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface ! Define the V grid FILE ( nid_V ) clhstnam = TRIM(cdfile_name)//".grid_V" IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & & nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdepv & "m", ipk, gdept_1d, nz_V, "down" ) ! ! Index of ocean points CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface ! ----------------------------------------------------------------- ! 2. Declare all the output fields as NETCDF variables ! ----------------------------------------------------------------- ! !!! nid_T : 3D !CALL histdef( nid_T, "votempis", "Insitu Temperature" , "C" , & ! ! & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! sshn & jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout ) CALL histdef( nid_T, "votempis", "Insitu Temperature" , "C" , & ! rinsitu_t & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) CALL histdef( nid_T, "maxu" , "Max Zonal Current" , "m/s" , & ! zwu & jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout ) CALL histdef( nid_T, "maxv" , "Max Meridional Current" , "m/s" , & ! zwv & jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout ) CALL histdef( nid_T, "maxz" , "Max Current Depth" , "m/s" , & ! zwz & jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout ) CALL histdef( nid_T, "sbt" , "Bottom Temperature" , "C" , & ! sbt & jpi, jpj, nh_T, 1 , 1, 1 , nz_T, 32, clop, zsto, zout ) CALL histend( nid_T, snc4chunks=snc4set ) ! !!! nid_U : 3D CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un & jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) CALL histend( nid_U, snc4chunks=snc4set ) ! !!! nid_V : 3D CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn & jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) CALL histend( nid_V, snc4chunks=snc4set ) ! ----------------------------------------------------------------- ! 3. Write the data ! ----------------------------------------------------------------- idex(1) = 1 CALL histwrite( nid_T, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex ) ! now temperature CALL histwrite( nid_T, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex ) ! now salinity CALL histwrite( nid_T, "sossheig", kt, sshn , jpi*jpj , idex ) ! sea surface height CALL theta2t CALL histwrite( nid_T, "votempis", kt, rinsitu_t(:,:,:) , jpi*jpj*jpk, idex ) ! now insitu-temperature CALL calc_max_cur(zwu,zwv,zwz,zmdi) CALL lbc_lnk( zwu, 'T', 1. ) CALL lbc_lnk( zwv, 'T', 1. ) CALL lbc_lnk( zwz, 'T', 1. ) CALL histwrite( nid_T, "maxu" , kt, zwu , jpi*jpj , idex ) ! max u-current CALL histwrite( nid_T, "maxv" , kt, zwv , jpi*jpj , idex ) ! max v-current CALL histwrite( nid_T, "maxz" , kt, zwz , jpi*jpj , idex ) ! max current depth DO jj = 1, jpj DO ji = 1, jpi jkbot = mbkt(ji,jj) z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem) END DO END DO CALL histwrite( nid_T, "sbt" , kt, z2d , jpi*jpj , idex ) ! sbt ! U file CALL histwrite( nid_U, "vozocrtx", kt, un , jpi*jpj*jpk, idex ) ! now i-velocity ! V file CALL histwrite( nid_V, "vomecrty", kt, vn , jpi*jpj*jpk, idex ) ! now j-velocity ! ----------------------------------------------------------------- ! 4. Close the files ! ----------------------------------------------------------------- CALL histclo( nid_T ) CALL histclo( nid_U ) CALL histclo( nid_V ) END SUBROUTINE dia_diaopfoam_zero END MODULE diaopfoam