Back to home page

MITgcm

 
 

    


File indexing completed on 2024-03-02 06:10:20 UTC

view on githubraw file Latest commit 5cf43646 on 2024-03-01 18:50:49 UTC
7bfe6112e8 Jean*0001 #include "CTRL_OPTIONS.h"
5d5c0b0d52 Patr*0002 
                0003       subroutine ctrl_set_pack_xz(
5cf4364659 Mart*0004      &     cunit, ivar, fname, masktype,weighttype,
5d5c0b0d52 Patr*0005      &     weightfld, lxxadxx, mythid)
                0006 
                0007 c     ==================================================================
                0008 c     SUBROUTINE ctrl_set_pack_xz
                0009 c     ==================================================================
                0010 c
                0011 c     o Compress the control vector such that only ocean points are
                0012 c       written to file.
                0013 c
7109a141b2 Patr*0014 c     o Open boundary packing finalized :
                0015 c          gebbie@mit.edu, 18-Mar-2003
                0016 c
                0017 c     changed: heimbach@mit.edu 17-Jun-2003
2146dab1aa Jean*0018 c              merged changes from Armin to replace write of
7109a141b2 Patr*0019 c              nr * globfld2d by 1 * globfld3d
                0020 c              (ad hoc fix to speed up global I/O)
                0021 c
5d5c0b0d52 Patr*0022 c     ==================================================================
                0023 
                0024       implicit none
                0025 
                0026 c     == global variables ==
                0027 
                0028 #include "EEPARAMS.h"
                0029 #include "SIZE.h"
                0030 #include "PARAMS.h"
                0031 #include "GRID.h"
                0032 
5cf4364659 Mart*0033 #include "CTRL_SIZE.h"
4d72283393 Mart*0034 #include "CTRL.h"
e612621177 Gael*0035 #include "CTRL_OBCS.h"
65754df434 Mart*0036 #include "OPTIMCYCLE.h"
5d5c0b0d52 Patr*0037 
                0038 c     == routine arguments ==
                0039 
                0040       integer cunit
5cf4364659 Mart*0041       integer ivar
de57a2ec4b Mart*0042       character*(MAX_LEN_FNAM) fname
5d5c0b0d52 Patr*0043       character*(  9) masktype
5cf4364659 Mart*0044       character*(  *) weighttype
de57a2ec4b Mart*0045       _RL     weightfld( Nr,nobcs )
5d5c0b0d52 Patr*0046       logical lxxadxx
                0047       integer mythid
                0048 
dac57cef35 Patr*0049 #ifndef EXCLUDE_CTRL_PACK
5d5c0b0d52 Patr*0050 c     == local variables ==
                0051 
                0052       integer bi,bj
                0053       integer ip,jp
                0054       integer i,j,k
abacc7d9db Patr*0055       integer ii,jj,kk
7109a141b2 Patr*0056       integer irec,iobcs,nrec_nl
5d5c0b0d52 Patr*0057       integer itlo,ithi
                0058       integer jtlo,jthi
                0059       integer jmin,jmax
                0060       integer imin,imax
                0061 
                0062       integer cbuffindex
f5224d0b03 Patr*0063       integer reclen, irectrue
                0064       integer cunit2, cunit3
de57a2ec4b Mart*0065       character*(MAX_LEN_FNAM) cfile2, cfile3
                0066 
                0067       real*4  cbuff      ( sNx*nSx*nPx*nSy*nPy )
                0068       real*4  globfldtmp2( sNx,nSx,nPx,nSy,nPy )
                0069       real*4  globfldtmp3( sNx,nSx,nPx,nSy,nPy )
                0070       _RL     globfldxz  ( sNx,nSx,nPx,nSy,nPy,Nr )
                0071       _RL     globfld3d  ( sNx,nSx,nPx,sNy,nSy,nPy,Nr )
                0072       _RL     globmskxz  ( sNx,nSx,nPx,nSy,nPy,Nr,nobcs )
7109a141b2 Patr*0073 #ifdef CTRL_PACK_PRECISE
5b7a03205a Mart*0074       integer il
de57a2ec4b Mart*0075       character*(MAX_LEN_FNAM) weightname
                0076       _RL     weightfldxz( sNx,nSx,nPx,nSy,nPy,Nr,nobcs )
5d5c0b0d52 Patr*0077 
                0078 c     == external ==
                0079 
                0080       integer  ilnblnk
                0081       external ilnblnk
5b7a03205a Mart*0082 #endif
5d5c0b0d52 Patr*0083 
                0084 c     == end of interface ==
                0085 
                0086       jtlo = 1
de57a2ec4b Mart*0087       jthi = nSy
5d5c0b0d52 Patr*0088       itlo = 1
de57a2ec4b Mart*0089       ithi = nSx
5d5c0b0d52 Patr*0090       jmin = 1
de57a2ec4b Mart*0091       jmax = sNy
5d5c0b0d52 Patr*0092       imin = 1
de57a2ec4b Mart*0093       imax = sNx
5d5c0b0d52 Patr*0094 
                0095 c     Initialise temporary file
de57a2ec4b Mart*0096       do k = 1,Nr
5d5c0b0d52 Patr*0097          do jp = 1,nPy
                0098             do bj = jtlo,jthi
                0099                do ip = 1,nPx
                0100                   do bi = itlo,ithi
                0101                      do i = imin,imax
f5224d0b03 Patr*0102                         globfldxz  (i,bi,ip,bj,jp,k) = 0. _d 0
                0103                         globfldtmp2(i,bi,ip,bj,jp)   = 0.
                0104                         globfldtmp3(i,bi,ip,bj,jp)   = 0.
7109a141b2 Patr*0105                         do iobcs=1,nobcs
                0106                            globmskxz(i,bi,ip,bj,jp,k,iobcs) = 0. _d 0
                0107                         enddo
                0108                      enddo
                0109                   enddo
                0110                enddo
                0111             enddo
                0112          enddo
                0113       enddo
                0114 c     Initialise temporary file
de57a2ec4b Mart*0115       do k = 1,Nr
7109a141b2 Patr*0116          do jp = 1,nPy
                0117             do bj = jtlo,jthi
                0118                do j = jmin,jmax
                0119                   do ip = 1,nPx
                0120                      do bi = itlo,ithi
                0121                         do i = imin,imax
                0122                            globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0
                0123                         enddo
5d5c0b0d52 Patr*0124                      enddo
                0125                   enddo
                0126                enddo
                0127             enddo
                0128          enddo
                0129       enddo
                0130 
                0131 c--   Only the master thread will do I/O.
                0132       _BEGIN_MASTER( mythid )
                0133 
f5224d0b03 Patr*0134       if ( doPackDiag ) then
                0135          if ( lxxadxx ) then
de57a2ec4b Mart*0136             write(cfile2,'(a,I2.2,a,I4.4,a)')
f5224d0b03 Patr*0137      &           'diag_pack_nonout_ctrl_',
5cf4364659 Mart*0138      &           ivar, '_', optimcycle, '.bin'
de57a2ec4b Mart*0139             write(cfile3,'(a,I2.2,a,I4.4,a)')
f5224d0b03 Patr*0140      &           'diag_pack_dimout_ctrl_',
5cf4364659 Mart*0141      &           ivar, '_', optimcycle, '.bin'
f5224d0b03 Patr*0142          else
de57a2ec4b Mart*0143             write(cfile2,'(a,I2.2,a,I4.4,a)')
f5224d0b03 Patr*0144      &           'diag_pack_nonout_grad_',
5cf4364659 Mart*0145      &           ivar, '_', optimcycle, '.bin'
de57a2ec4b Mart*0146             write(cfile3,'(a,I2.2,a,I4.4,a)')
f5224d0b03 Patr*0147      &           'diag_pack_dimout_grad_',
5cf4364659 Mart*0148      &           ivar, '_', optimcycle, '.bin'
f5224d0b03 Patr*0149          endif
                0150 
de57a2ec4b Mart*0151          reclen = sNx*nSx*nPx*nSy*nPy*4
f5224d0b03 Patr*0152          call mdsfindunit( cunit2, mythid )
                0153          open( cunit2, file=cfile2, status='unknown',
                0154      &        access='direct', recl=reclen )
                0155          call mdsfindunit( cunit3, mythid )
                0156          open( cunit3, file=cfile3, status='unknown',
                0157      &        access='direct', recl=reclen )
                0158       endif
                0159 
7109a141b2 Patr*0160       do iobcs = 1, nobcs
                0161          call MDSREADFIELD_XZ_GL(
                0162      &        masktype, ctrlprec, 'RL',
                0163      &        Nr, globmskxz(1,1,1,1,1,1,iobcs), iobcs, mythid)
                0164 #ifdef CTRL_PACK_PRECISE
                0165          il=ilnblnk( weighttype)
de57a2ec4b Mart*0166          write(weightname,'(a)') weighttype(1:il)
7109a141b2 Patr*0167          call MDSREADFIELD_XZ_GL(
                0168      &        weightname, ctrlprec, 'RL',
                0169      &        Nr, weightfldxz(1,1,1,1,1,1,iobcs), iobcs, mythid)
                0170 #endif
                0171       enddo
                0172 
b94917e7bd Mart*0173       if ( useSingleCPUio ) then
                0174 C     MDSREADFIELD_XZ_GL does not know about useSingleCPUio, so the faster
                0175 C     method that works for .not.useSingleCPUio cannot be used
                0176        nrec_nl = 0
                0177       else
5cf4364659 Mart*0178        nrec_nl = int(ncvarrecs(ivar)/Ny)
b94917e7bd Mart*0179       endif
7109a141b2 Patr*0180       do irec = 1, nrec_nl
                0181          call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
de57a2ec4b Mart*0182      &        Nr, globfld3d, irec, mythid)
                0183          do j=1,sNy
                0184             iobcs= mod((irec-1)*sNy+j-1,nobcs)+1
7109a141b2 Patr*0185 
5cf4364659 Mart*0186             write(cunit) ncvarindex(ivar)
7109a141b2 Patr*0187             write(cunit) 1
                0188             write(cunit) 1
de57a2ec4b Mart*0189             do k = 1,Nr
                0190              irectrue = (irec-1)*nobcs*Nr + (iobcs-1)*Nr + k
7109a141b2 Patr*0191              cbuffindex = 0
                0192              do jp = 1,nPy
                0193               do bj = jtlo,jthi
                0194                do ip = 1,nPx
                0195                 do bi = itlo,ithi
                0196                  do i = imin,imax
de57a2ec4b Mart*0197                   jj=mod((j-1)*Nr+k-1,sNy)+1
                0198                   kk=int((j-1)*Nr+K-1)/sNy+1
7109a141b2 Patr*0199                   if (globmskxz(i,bi,ip,bj,jp,k,iobcs)  .ne. 0. ) then
                0200                      cbuffindex = cbuffindex + 1
f5224d0b03 Patr*0201 cph(
                0202                      globfldtmp3(i,bi,ip,bj,jp) =
                0203      &                    globfld3d(i,bi,ip,jj,bj,jp,kk)
                0204 cph)
7109a141b2 Patr*0205 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                0206                      if (lxxadxx) then
951926fb9b Jean*0207                         cbuff(cbuffindex) =
abacc7d9db Patr*0208      &                       globfld3d(i,bi,ip,jj,bj,jp,kk) *
7109a141b2 Patr*0209 # ifdef CTRL_PACK_PRECISE
                0210      &                       sqrt(weightfldxz(i,bi,ip,bj,jp,k,iobcs))
                0211 # else
                0212      &                       sqrt(weightfld(k,iobcs))
                0213 # endif
                0214                      else
951926fb9b Jean*0215                         cbuff(cbuffindex) =
abacc7d9db Patr*0216      &                       globfld3d(i,bi,ip,jj,bj,jp,kk) /
7109a141b2 Patr*0217 # ifdef CTRL_PACK_PRECISE
                0218      &                       sqrt(weightfldxz(i,bi,ip,bj,jp,k,iobcs))
                0219 # else
                0220      &                       sqrt(weightfld(k,iobcs))
                0221 # endif
                0222                      endif
f5224d0b03 Patr*0223 cph(
                0224                      globfldtmp2(i,bi,ip,bj,jp) = cbuff(cbuffindex)
                0225 cph)
7109a141b2 Patr*0226 #else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
abacc7d9db Patr*0227                      cbuff(cbuffindex) = globfld3d(i,bi,ip,jj,bj,jp,kk)
7109a141b2 Patr*0228 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
                0229                   endif
                0230                  enddo
                0231                 enddo
                0232                enddo
                0233               enddo
                0234              enddo
                0235 c           --> check cbuffindex.
                0236              if ( cbuffindex .gt. 0) then
                0237                 write(cunit) cbuffindex
                0238                 write(cunit) k
                0239                 write(cunit) (cbuff(ii), ii=1,cbuffindex)
                0240              endif
f5224d0b03 Patr*0241 c
                0242             if ( doPackDiag ) then
                0243                write(cunit2,rec=irectrue) globfldtmp2
                0244                write(cunit3,rec=irectrue) globfldtmp3
                0245             endif
                0246 c
7109a141b2 Patr*0247 c     -- end of k loop --
                0248             enddo
                0249 c     -- end of j loop --
                0250          enddo
                0251 c     -- end of irec loop --
                0252       enddo
                0253 
5cf4364659 Mart*0254       do irec = nrec_nl*ny+1, ncvarrecs(ivar)
5b7a03205a Mart*0255 c     Need to solve for what iobcs would have been.
7109a141b2 Patr*0256          iobcs= mod(irec-1,nobcs)+1
5d5c0b0d52 Patr*0257 
                0258          call MDSREADFIELD_XZ_GL( fname, ctrlprec, 'RL',
de57a2ec4b Mart*0259      &        Nr, globfldxz, irec, mythid)
5d5c0b0d52 Patr*0260 
5cf4364659 Mart*0261          write(cunit) ncvarindex(ivar)
5d5c0b0d52 Patr*0262          write(cunit) 1
                0263          write(cunit) 1
de57a2ec4b Mart*0264          do k = 1,Nr
                0265             irectrue = (irec-1)*nobcs*Nr + (iobcs-1)*Nr + k
5d5c0b0d52 Patr*0266             cbuffindex = 0
                0267             do jp = 1,nPy
                0268              do bj = jtlo,jthi
                0269               do ip = 1,nPx
                0270                do bi = itlo,ithi
                0271                 do i = imin,imax
7109a141b2 Patr*0272                  if (globmskxz(i,bi,ip,bj,jp,k,iobcs)  .ne. 0. ) then
5d5c0b0d52 Patr*0273                      cbuffindex = cbuffindex + 1
f5224d0b03 Patr*0274 cph(
                0275                      globfldtmp3(i,bi,ip,bj,jp) =
                0276      &                    globfldxz(i,bi,ip,bj,jp,k)
                0277 cph)
5d5c0b0d52 Patr*0278 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                0279                      if (lxxadxx) then
951926fb9b Jean*0280                         cbuff(cbuffindex) =
5d5c0b0d52 Patr*0281      &                       globfldxz(i,bi,ip,bj,jp,k) *
7109a141b2 Patr*0282 # ifdef CTRL_PACK_PRECISE
                0283      &                       sqrt(weightfldxz(i,bi,ip,bj,jp,k,iobcs))
                0284 # else
5d5c0b0d52 Patr*0285      &                       sqrt(weightfld(k,iobcs))
7109a141b2 Patr*0286 # endif
5d5c0b0d52 Patr*0287                      else
951926fb9b Jean*0288                         cbuff(cbuffindex) =
5d5c0b0d52 Patr*0289      &                       globfldxz(i,bi,ip,bj,jp,k) /
7109a141b2 Patr*0290 # ifdef CTRL_PACK_PRECISE
                0291      &                       sqrt(weightfldxz(i,bi,ip,bj,jp,k,iobcs))
                0292 # else
5d5c0b0d52 Patr*0293      &                       sqrt(weightfld(k,iobcs))
7109a141b2 Patr*0294 # endif
5d5c0b0d52 Patr*0295                      endif
f5224d0b03 Patr*0296 cph(
                0297                      globfldtmp2(i,bi,ip,bj,jp) = cbuff(cbuffindex)
                0298 cph)
                0299 #else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
5d5c0b0d52 Patr*0300                      cbuff(cbuffindex) = globfldxz(i,bi,ip,bj,jp,k)
f5224d0b03 Patr*0301 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
5d5c0b0d52 Patr*0302                  endif
                0303                 enddo
                0304                enddo
                0305               enddo
                0306              enddo
                0307             enddo
                0308 c           --> check cbuffindex.
                0309             if ( cbuffindex .gt. 0) then
                0310                write(cunit) cbuffindex
                0311                write(cunit) k
                0312                write(cunit) (cbuff(ii), ii=1,cbuffindex)
                0313             endif
                0314 c
f5224d0b03 Patr*0315             if ( doPackDiag ) then
                0316                write(cunit2,rec=irectrue) globfldtmp2
                0317                write(cunit3,rec=irectrue) globfldtmp3
                0318             endif
                0319 c
7109a141b2 Patr*0320 c     -- end of k loop --
                0321          enddo
5d5c0b0d52 Patr*0322 c     -- end of irec loop --
                0323       enddo
                0324 
                0325       _END_MASTER( mythid )
                0326 
dac57cef35 Patr*0327 #endif
                0328 
5d5c0b0d52 Patr*0329       return
                0330       end