MODULE domwri !!====================================================================== !! *** MODULE domwri *** !! Ocean initialization : write the ocean domain mesh ask file(s) !!====================================================================== !!---------------------------------------------------------------------- !! dom_wri : create and write mesh and mask file(s) !! nmsh = 1 : mesh_mask file !! = 2 : mesh and mask file !! = 3 : mesh_hgr, mesh_zgr and mask !!---------------------------------------------------------------------- !! * Modules used USE dom_oce ! ocean space and time domain USE in_out_manager USE iom USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE lib_mpp IMPLICIT NONE PRIVATE PUBLIC dom_wri ! routine called by inidom.F90 !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Id: domwri.F90 1590 2009-08-06 10:18:30Z smasson $ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dom_wri !!---------------------------------------------------------------------- !! *** ROUTINE dom_wri *** !! !! ** Purpose : Create the NetCDF file(s) which contain(s) all the !! ocean domain informations (mesh and mask arrays). This (these) !! file(s) is (are) used for visualisation (SAXO software) and !! diagnostic computation. !! !! ** Method : Write in a file all the arrays generated in routines !! domhgr, domzgr, and dommsk. Note: the file contain depends on !! the vertical coord. used (z-coord, partial steps, s-coord) !! nmsh = 1 : 'mesh_mask.nc' file !! = 2 : 'mesh.nc' and mask.nc' files !! = 3 : 'mesh_hgr.nc', 'mesh_zgr.nc' and !! 'mask.nc' files !! For huge size domain, use option 2 or 3 depending on your !! vertical coordinate. !! !! ** output file : !! meshmask.nc : domain size, horizontal grid-point position, !! masks, depth and vertical scale factors !! !! History : !! ! 97-02 (G. Madec) Original code !! ! 99-11 (M. Imbard) NetCDF FORMAT with IOIPSL !! 9.0 ! 02-08 (G. Madec) F90 and several file !!---------------------------------------------------------------------- INTEGER :: inum0 ! temprary units for 'mesh_mask.nc' file INTEGER :: inum1 ! temprary units for 'mesh.nc' file INTEGER :: inum2 ! temprary units for 'mask.nc' file INTEGER :: inum3 ! temprary units for 'mesh_hgr.nc' file INTEGER :: inum4 ! temprary units for 'mesh_zgr.nc' file INTEGER :: ji, jj, jk, ik REAL(wp), DIMENSION(jpi,jpj) :: zprt ! temporary array for bathymetry REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepu ! 3D depth of U point REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepv ! 3D depth of V point CHARACTER(len=21) :: clnam0 ! filename (mesh and mask informations) CHARACTER(len=21) :: clnam1 ! filename (mesh informations) CHARACTER(len=21) :: clnam2 ! filename (mask informations) CHARACTER(len=21) :: clnam3 ! filename (horizontal mesh informations) CHARACTER(len=21) :: clnam4 ! filename (vertical mesh informations) !!---------------------------------------------------------------------- IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dom_wri : create NetCDF mesh and mask information file(s)' IF(lwp) WRITE(numout,*) '~~~~~~~' clnam0 = 'mesh_mask' ! filename (mesh and mask informations) clnam1 = 'mesh' ! filename (mesh informations) clnam2 = 'mask' ! filename (mask informations) clnam3 = 'mesh_hgr' ! filename (horizontal mesh informations) clnam4 = 'mesh_zgr' ! filename (vertical mesh informations) SELECT CASE ( MOD(nmsh, 3) ) ! ! ============================ CASE ( 1 ) ! create 'mesh_mask.nc' file ! ! ============================ CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib ) inum2 = inum0 ! put all the informations inum3 = inum0 ! in unit inum0 inum4 = inum0 ! ! ============================ CASE ( 2 ) ! create 'mesh.nc' and ! ! 'mask.nc' files ! ! ============================ CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib ) CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib ) inum3 = inum1 ! put mesh informations inum4 = inum1 ! in unit inum1 ! ! ============================ CASE ( 0 ) ! create 'mesh_hgr.nc' ! ! 'mesh_zgr.nc' and ! ! 'mask.nc' files ! ! ============================ CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib ) CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib ) CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib ) END SELECT ! ! masks (inum2) CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 ) ! ! land-sea mask CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 ) CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 ) CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 ) zprt = tmask(:,:,1) * dom_uniq('T') ! ! unique point mask CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) zprt = umask(:,:,1) * dom_uniq('U') CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) zprt = vmask(:,:,1) * dom_uniq('V') CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) zprt = fmask(:,:,1) * dom_uniq('F') CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) ! ! horizontal mesh (inum3) CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r4 ) ! ! latitude CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r4 ) ! ! longitude CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 ) CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 ) CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 ) CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 ) ! ! e2 scale factors CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 ) CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 ) CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 ) CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 ) ! ! coriolis factor ! note that mbathy has been modified in dommsk or in solver. ! it is the number of non-zero "w" levels in the water, and the minimum ! value (on land) is 2. We define zprt as the number of "T" points in the ocean ! at any location, and zero on land. ! zprt = tmask(:,:,1)*(mbathy-1) CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 ) #if ! defined key_zco IF( ln_sco ) THEN ! s-coordinate CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) ! ! depth CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt ) ! ! scaling coef. CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w ) CALL iom_rstput( 0, 0, inum4, 'esigt', esigt ) CALL iom_rstput( 0, 0, inum4, 'esigw', esigw ) CALL iom_rstput( 0, 0, inum4, 'e3t', e3t ) ! ! scale factors CALL iom_rstput( 0, 0, inum4, 'e3u', e3u ) CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 ) ! ! stretched system CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) ENDIF IF( ln_zps ) THEN ! z-coordinate - partial steps IF( nmsh <= 6 ) THEN ! ! 3D vertical scale factors CALL iom_rstput( 0, 0, inum4, 'e3t', e3t ) CALL iom_rstput( 0, 0, inum4, 'e3u', e3u ) CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) ELSE ! ! 2D bottom scale factors DO jj = 1,jpj ; DO ji = 1,jpi ik = NINT( zprt(ji,jj) ) ! take care that mbathy is not what you think it is here ! IF ( ik /= 0 ) THEN ; e3tp(ji,jj) = e3t(ji,jj,ik) ; e3wp(ji,jj) = e3w(ji,jj,ik) ELSE ; e3tp(ji,jj) = 0. ; e3wp(ji,jj) = 0. ENDIF END DO ; END DO CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp ) CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp ) END IF IF( nmsh <= 3 ) THEN ! ! 3D depth CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 ) DO jk = 1,jpk ; DO jj = 1, jpjm1 ; DO ji = 1, fs_jpim1 ! vector opt. zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji+1,jj ,jk) ) zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji ,jj+1,jk) ) END DO ; END DO ; END DO CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. ) CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 ) ELSE ! ! 2D bottom depth DO jj = 1,jpj ; DO ji = 1,jpi ik = NINT( zprt(ji,jj) ) ! take care that mbathy is not what you think it is here ! IF ( ik /= 0 ) THEN ; hdept(ji,jj) = gdept(ji,jj,ik) ; hdepw(ji,jj) = gdepw(ji,jj,ik+1) ELSE ; hdept(ji,jj) = 0. ; hdepw(ji,jj) = 0. ENDIF END DO ; END DO CALL iom_rstput( 0, 0, inum4, 'hdept' , hdept, ktype = jp_r4 ) CALL iom_rstput( 0, 0, inum4, 'hdepw' , hdepw, ktype = jp_r4 ) ENDIF CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 ) ! ! reference z-coord. CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 ) CALL iom_rstput( 0, 0, inum4, 'e3t_0' , e3t_0 ) CALL iom_rstput( 0, 0, inum4, 'e3w_0' , e3w_0 ) ENDIF #endif IF( ln_zco ) THEN ! ! z-coordinate - full steps CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 ) ! ! depth CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 ) CALL iom_rstput( 0, 0, inum4, 'e3t_0' , e3t_0 ) ! ! scale factors CALL iom_rstput( 0, 0, inum4, 'e3w_0' , e3w_0 ) ENDIF ! ! ============================ ! ! close the files ! ! ============================ SELECT CASE ( nmsh ) CASE ( 1 ) CALL iom_close( inum0 ) CASE ( 2 ) CALL iom_close( inum1 ) CALL iom_close( inum2 ) CASE ( 3 ) CALL iom_close( inum2 ) CALL iom_close( inum3 ) CALL iom_close( inum4 ) END SELECT END SUBROUTINE dom_wri FUNCTION dom_uniq( cdgrd ) RESULT( puniq ) !!---------------------------------------------------------------------- !! *** ROUTINE dom_uniq *** !! !! ** Purpose : identify unique point of a grid (TUVF) !! !! ** Method : 1) aplly lbc_lnk on an array with different values for each element !! 2) check which elements have been changed !! !!---------------------------------------------------------------------- CHARACTER(len=1) , INTENT(in ) :: cdgrd ! REAL(wp), DIMENSION(jpi,jpj) :: puniq ! REAL(wp), DIMENSION(jpi,jpj ) :: ztstref ! array with different values for each element REAL(wp) :: zshift ! shift value link to the process number LOGICAL , DIMENSION(jpi,jpj,1) :: lldbl ! is the point unique or not? INTEGER :: ji ! dummy loop indices !!---------------------------------------------------------------------- ! build an array with different values for each element ! in mpp: make sure that these values are different even between process ! -> apply a shift value according to the process number zshift = jpi * jpj * ( narea - 1 ) ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji, wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) ) puniq(:,:) = ztstref(:,:) ! default definition CALL lbc_lnk( puniq, cdgrd, 1. ) ! apply boundary conditions lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed puniq(:,:) = 1. ! default definition ! fill only the inner part of the cpu with llbl converted into real puniq(nldi:nlei,nldj:nlej) = REAL(COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ), wp) END FUNCTION dom_uniq !!====================================================================== END MODULE domwri