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 12536 for branches – NEMO

Changeset 12536 for branches


Ignore:
Timestamp:
2020-03-11T14:50:30+01:00 (4 years ago)
Author:
petesykes
Message:

added diaopfoam T0 diagnostics

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_PS44/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90

    r10390 r12536  
    1515   USE diurnal_bulk 
    1616   USE cool_skin 
     17   USE ioipsl 
    1718 
    1819 
     
    2122 
    2223   LOGICAL , PUBLIC ::   ln_diaopfoam     !: Diaopfoam output  
     24   LOGICAL , PUBLIC ::   ln_diaopfoam_Tzero     !: Diaopfoam first time step output 
    2325   PUBLIC   dia_diaopfoam_init            ! routine called by nemogcm.F90 
    2426   PUBLIC   dia_diaopfoam                 ! routine called by diawri.F90 
     
    4951      INTEGER            ::   ierror   ! local integer 
    5052      !! 
    51       NAMELIST/nam_diadiaop/ ln_diaopfoam 
     53      NAMELIST/nam_diadiaop/ ln_diaopfoam,ln_diaopfoam_Tzero 
    5254      !! 
    5355      !!---------------------------------------------------------------------- 
    5456      ! 
    5557      ln_diaopfoam = .false.         ! default value for diaopfoam stream 
     58      ln_diaopfoam_Tzero = .false.         ! default value for diaopfoam Tzero stream 
    5659      REWIND ( numnam_ref )              ! Read Namelist nam_diadiaop in reference namelist : 3D hourly diagnostics 
    5760      READ   ( numnam_ref, nam_diadiaop, IOSTAT=ios, ERR= 901 ) 
     
    6871         WRITE(numout,*) '~~~~~~~~~~~~' 
    6972         WRITE(numout,*) '   Namelist nam_diadiaop : set diaopfoam outputs ' 
    70          WRITE(numout,*) '      Switch for diaopfoam diagnostics (T) or not (F)  ln_diaopfoam  = ', ln_diaopfoam 
     73         WRITE(numout,*) '      Switch for diaopfoam diagnostics (T) or not (F)                 ln_diaopfoam        = ', ln_diaopfoam 
     74         WRITE(numout,*) '      Switch for diaopfoam first timestep diagnostics (T) or not (F)  ln_diaopfoam_Tzero  = ', ln_diaopfoam_Tzero 
    7175      ENDIF 
    7276    END SUBROUTINE dia_diaopfoam_init 
    7377 
    74     SUBROUTINE dia_diaopfoam 
     78    SUBROUTINE dia_diaopfoam( kt ) 
    7579      !!---------------------------------------------------------------------- 
    7680      !!                 ***  ROUTINE dia_diaopfoam  *** 
     
    8286      !!          
    8387      !!-------------------------------------------------------------------- 
     88      IMPLICIT NONE 
     89 
     90      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     91 
    8492      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace 
    8593      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
     
    94102 
    95103      zmdi=1.e+20                               !  missing data indicator for masking 
     104 
    96105      ! Diaopfoam stream if needed 
    97106      IF (ln_diaopfoam) THEN 
     107         IF ( kt .eq. nn_it000 .AND. ln_diaopfoam_Tzero ) THEN 
     108            IF(lwp) WRITE(numout,*) 'diaopfoam: writing T0 at kt = ', kt 
     109            CALL dia_diaopfoam_zero( 'Tzero', kt ) 
     110         ENDIF 
     111          
    98112         CALL theta2t 
    99113         CALL iom_put( "insitut_op" , rinsitu_t(:,:,:)                      )    ! insitu temperature 
     
    101115         CALL iom_put( "soce_op"   , tsn(:,:,:,jp_sal)                     )    ! salinity 
    102116         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  
     117            CALL iom_put( "sst_wl_op"   , x_dsst                             )    ! warm layer  
     118            CALL iom_put( "sst_cs_op"   , x_csdsst                           )    ! cool skin  
    105119         ENDIF 
    106120         zw2d(:,:)=sshn(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     
    175189   END SUBROUTINE calc_max_cur 
    176190 
     191   SUBROUTINE dia_diaopfoam_zero( cdfile_name, kt ) 
     192      !!--------------------------------------------------------------------- 
     193      !!                 ***  ROUTINE dia_diaopfoam_zero  *** 
     194      !!         
     195      !! ** Purpose :   create a NetCDF file named cdfile_name which contains  
     196      !!      the instantaneous ocean state at the first tiome step. 
     197      !! 
     198      !! ** Method  :   NetCDF files using ioipsl 
     199      !!---------------------------------------------------------------------- 
     200      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created 
     201      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index 
     202      !!  
     203      CHARACTER (len=32) :: clhstnam 
     204      CHARACTER (len=40) :: clop 
     205      INTEGER  ::   iimi, iima, ipk, ijmi, ijma          ! local integers 
     206      INTEGER  ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file 
     207      INTEGER  ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file 
     208      INTEGER  ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
     209      INTEGER  ::   id_i , nz_i, nh_i        
     210      INTEGER, DIMENSION(1) ::   idex                    ! local workspace 
     211      INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
     212      INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
     213      INTEGER :: ierr 
     214      INTEGER :: jkbot, jj, ji 
     215      REAL(wp) :: zsto, zout, zmax, zjulian, zdt 
     216      REAL(wp) :: zmdi 
     217      REAL(wp), POINTER, DIMENSION(:,:)   :: zwu 
     218      REAL(wp), POINTER, DIMENSION(:,:)   :: zwv 
     219      REAL(wp), POINTER, DIMENSION(:,:)   :: zwz 
     220      REAL(wp), DIMENSION(jpi,jpj) :: zw2d           ! 2D workspace 
     221      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
     222      REAL(wp), DIMENSION(jpi,jpj) :: z2d 
     223 
     224      !!---------------------------------------------------------------------- 
     225 
     226      ! ----------------------------------------------------------------- 
     227      ! 0. Allocations 
     228      ! ----------------------------------------------------------------- 
     229      ierr = 0 
     230      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
     231         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
     232         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr ) 
     233      IF( lk_mpp )   CALL mpp_sum( ierr ) 
     234         IF( ierr /= 0 ) THEN 
     235            CALL ctl_stop('dia_diaopfoam_zero: failed to allocate arrays') 
     236            RETURN 
     237      ENDIF 
     238      CALL wrk_alloc( jpi , jpj      , zwu ) 
     239      CALL wrk_alloc( jpi , jpj      , zwv ) 
     240      CALL wrk_alloc( jpi , jpj      , zwz ) 
     241 
     242      zmdi=1.e+20                               !  missing data indicator for masking 
     243 
     244      ! ----------------------------------------------------------------- 
     245      ! 1. Define NETCDF files and fields at beginning of first time step 
     246      ! ----------------------------------------------------------------- 
     247 
     248      ! Define indices of the horizontal output zoom and vertical limit storage 
     249      iimi = 1      ;      iima = jpi 
     250      ijmi = 1      ;      ijma = jpj 
     251      ipk = jpk 
     252 
     253      ! Define frequency of output and means 
     254      zdt  = rdt 
     255      zsto = rdt 
     256      clop = "inst(x)"           ! no use of the mask value (require less cpu time) 
     257      zout = rdt 
     258      zmax = ( nitend - nit000 + 1 ) * zdt 
     259 
     260      ! Compute julian date from starting date of the run 
     261      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis  
     262      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment 
     263 
     264      ! Define the T grid FILE ( nid_T ) 
     265      clhstnam = TRIM(cdfile_name)//".grid_T" 
     266      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
     267      CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     268         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
     269         &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 
     270      CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept 
     271         &           "m", ipk, gdept_1d, nz_T, "down" ) 
     272      !                                                            ! Index of ocean points 
     273      CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume 
     274      CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface 
     275 
     276      ! Define the U grid FILE ( nid_U ) 
     277      clhstnam = TRIM(cdfile_name)//".grid_U" 
     278      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
     279      CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu 
     280         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
     281         &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 
     282      CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdepu 
     283         &           "m", ipk, gdept_1d, nz_U, "down" ) 
     284      !                                                            ! Index of ocean points 
     285      CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume 
     286      CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface 
     287 
     288      ! Define the V grid FILE ( nid_V ) 
     289      clhstnam = TRIM(cdfile_name)//".grid_V" 
     290      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam 
     291      CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv 
     292         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
     293         &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 
     294      CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdepv 
     295         &          "m", ipk, gdept_1d, nz_V, "down" ) 
     296      !                                                            ! Index of ocean points 
     297      CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume 
     298      CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface 
     299 
     300 
     301      ! ----------------------------------------------------------------- 
     302      ! 2. Declare all the output fields as NETCDF variables 
     303      ! ----------------------------------------------------------------- 
     304 
     305      !                                                                                     !!! nid_T : 3D 
     306      !CALL histdef( nid_T, "votempis", "Insitu Temperature" , "C"      ,       &  !  
     307      !   &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
     308      CALL histdef( nid_T, "votemper", "Temperature"         , "C"     ,       &  ! tn 
     309         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
     310      CALL histdef( nid_T, "vosaline", "Salinity"            , "PSU"   ,       & ! sn 
     311         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
     312      CALL histdef( nid_T, "sossheig", "Sea Surface Height"  , "m"     ,       &  ! sshn 
     313         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout ) 
     314      CALL histdef( nid_T, "votempis", "Insitu Temperature"  , "C"     ,       &  ! rinsitu_t 
     315         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
     316      CALL histdef( nid_T, "maxu"    , "Max Zonal Current"   , "m/s"     ,       &  ! zwu 
     317         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout ) 
     318      CALL histdef( nid_T, "maxv"    , "Max Meridional Current" , "m/s"     ,    &  ! zwv 
     319         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout ) 
     320      CALL histdef( nid_T, "maxz"    , "Max Current Depth"   , "m/s"     ,       &  ! zwz 
     321         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout ) 
     322      CALL histdef( nid_T, "sbt"       , "Bottom Temperature"  , "C"     ,       &  ! sbt 
     323         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout ) 
     324      CALL histend( nid_T, snc4chunks=snc4set ) 
     325 
     326      !                                                                                     !!! nid_U : 3D 
     327      CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un 
     328         &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 
     329      CALL histend( nid_U, snc4chunks=snc4set ) 
     330 
     331      !                                                                                     !!! nid_V : 3D 
     332      CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn 
     333         &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 
     334      CALL histend( nid_V, snc4chunks=snc4set ) 
     335 
     336 
     337      ! ----------------------------------------------------------------- 
     338      ! 3. Write the data 
     339      ! ----------------------------------------------------------------- 
     340 
     341      idex(1) = 1 
     342      CALL histwrite( nid_T, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature 
     343      CALL histwrite( nid_T, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity 
     344      CALL histwrite( nid_T, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height 
     345      CALL theta2t 
     346      CALL histwrite( nid_T, "votempis", kt, rinsitu_t(:,:,:) , jpi*jpj*jpk, idex )    ! now insitu-temperature 
     347      CALL calc_max_cur(zwu,zwv,zwz,zmdi) 
     348      CALL lbc_lnk( zwu, 'T', 1. ) 
     349      CALL lbc_lnk( zwv, 'T', 1. ) 
     350      CALL lbc_lnk( zwz, 'T', 1. ) 
     351      CALL histwrite( nid_T, "maxu"    , kt, zwu              , jpi*jpj    , idex )    ! max u-current 
     352      CALL histwrite( nid_T, "maxv"    , kt, zwv              , jpi*jpj    , idex )    ! max v-current 
     353      CALL histwrite( nid_T, "maxz"    , kt, zwz              , jpi*jpj    , idex )    ! max current depth 
     354      DO jj = 1, jpj 
     355         DO ji = 1, jpi 
     356            jkbot = mbkt(ji,jj) 
     357            z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem) 
     358         END DO 
     359      END DO 
     360      CALL histwrite( nid_T, "sbt"     , kt, z2d              , jpi*jpj    , idex )    ! sbt 
     361      ! U file 
     362      CALL histwrite( nid_U, "vozocrtx", kt, un              , jpi*jpj*jpk, idex )     ! now i-velocity 
     363      ! V file 
     364      CALL histwrite( nid_V, "vomecrty", kt, vn              , jpi*jpj*jpk, idex )     ! now j-velocity 
     365 
     366 
     367      ! ----------------------------------------------------------------- 
     368      ! 4. Close the files 
     369      ! ----------------------------------------------------------------- 
     370 
     371      CALL histclo( nid_T ) 
     372      CALL histclo( nid_U ) 
     373      CALL histclo( nid_V ) 
     374 
     375   END SUBROUTINE dia_diaopfoam_zero 
     376 
    177377END MODULE diaopfoam 
Note: See TracChangeset for help on using the changeset viewer.