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_yz(
5cf4364659 Mart*0004      &     cunit, ivar, fname, masktype, weighttype,
5d5c0b0d52 Patr*0005      &     weightfld, lxxadxx, mythid)
                0006 
                0007 c     ==================================================================
                0008 c     SUBROUTINE ctrl_set_pack_yz
                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 added :
                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
dff4940422 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
5b7a03205a Mart*0066 
de57a2ec4b Mart*0067       real*4  cbuff      ( nSx*nPx*sNy*nSy*nPy )
                0068       real*4  globfldtmp2( nSx,nPx,sNy,nSy,nPy )
                0069       real*4  globfldtmp3( nSx,nPx,sNy,nSy,nPy )
                0070       _RL     globfldyz  ( nSx,nPx,sNy,nSy,nPy,Nr )
                0071       _RL     globfld3d  ( sNx,nSx,nPx,sNy,nSy,nPy,Nr )
                0072       _RL     globmskyz  ( nSx,nPx,sNy,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     weightfldyz( nSx,nPx,sNy,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 j = jmin,jmax
                0100                   do ip = 1,nPx
                0101                      do bi = itlo,ithi
f5224d0b03 Patr*0102                         globfldyz  (bi,ip,j,bj,jp,k) = 0. _d 0
                0103                         globfldtmp2(bi,ip,j,bj,jp)   = 0.
                0104                         globfldtmp3(bi,ip,j,bj,jp)   = 0.
7109a141b2 Patr*0105                         do iobcs=1,nobcs
                0106                            globmskyz(bi,ip,j,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 = nSx*nPx*sNy*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
951926fb9b Jean*0161          call MDSREADFIELD_YZ_GL(
7109a141b2 Patr*0162      &        masktype, ctrlprec, 'RL',
                0163      &        Nr, globmskyz(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_YZ_GL(
                0168      &       weightname, ctrlprec, 'RL',
                0169      &       Nr, weightfldyz(1,1,1,1,1,1,iobcs), iobcs, mythid)
                0170 #endif
                0171       enddo
5d5c0b0d52 Patr*0172 
b94917e7bd Mart*0173       if ( useSingleCPUio ) then
                0174 C     MDSREADFIELD_YZ_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)/Nx)
b94917e7bd Mart*0179       endif
7109a141b2 Patr*0180       do irec = 1, nrec_nl
5b7a03205a Mart*0181 c     Need to solve for what iobcs would have been.
5d5c0b0d52 Patr*0182 
7109a141b2 Patr*0183          call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
de57a2ec4b Mart*0184      &        Nr, globfld3D, irec, mythid)
7109a141b2 Patr*0185 
de57a2ec4b Mart*0186          do i=1,sNx
                0187             iobcs= mod((irec-1)*sNx+i-1,nobcs)+1
7109a141b2 Patr*0188 
5cf4364659 Mart*0189             write(cunit) ncvarindex(ivar)
7109a141b2 Patr*0190             write(cunit) 1
                0191             write(cunit) 1
de57a2ec4b Mart*0192             do k = 1,Nr
                0193              irectrue = (irec-1)*nobcs*Nr + (iobcs-1)*Nr + k
7109a141b2 Patr*0194              cbuffindex = 0
                0195              do jp = 1,nPy
                0196               do bj = jtlo,jthi
                0197                do ip = 1,nPx
                0198                 do bi = itlo,ithi
                0199                  do j = jmin,jmax
de57a2ec4b Mart*0200                   ii=mod ( (i-1)*Nr*sNy+(k-1)*sNy+j-1      , sNx ) + 1
                0201                   jj=mod( ((i-1)*Nr*sNy+(k-1)*sNy+j-1)/sNx , sNy ) + 1
                0202                   kk=int((i-1)*Nr*sNy+(k-1)*sNy+j-1)/(sNx*sNy) + 1
7109a141b2 Patr*0203                   if (globmskyz(bi,ip,j,bj,jp,k,iobcs)  .ne. 0. ) then
                0204                      cbuffindex = cbuffindex + 1
f5224d0b03 Patr*0205 cph(
                0206                      globfldtmp3(bi,ip,j,bj,jp) =
                0207      &                    globfld3d(ii,bi,ip,jj,bj,jp,kk)
                0208 cph)
7109a141b2 Patr*0209 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                0210                      if (lxxadxx) then
951926fb9b Jean*0211                         cbuff(cbuffindex) =
dff4940422 Patr*0212      &                       globfld3d(ii,bi,ip,jj,bj,jp,kk) *
7109a141b2 Patr*0213 #ifdef CTRL_PACK_PRECISE
                0214      &                       sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
                0215 #else
                0216      &                       sqrt(weightfld(k,iobcs))
                0217 #endif
                0218                      else
951926fb9b Jean*0219                         cbuff(cbuffindex) =
dff4940422 Patr*0220      &                       globfld3d(ii,bi,ip,jj,bj,jp,kk) /
7109a141b2 Patr*0221 #ifdef CTRL_PACK_PRECISE
                0222      &                       sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
                0223 #else
                0224      &                       sqrt(weightfld(k,iobcs))
                0225 #endif
                0226                      endif
f5224d0b03 Patr*0227 cph(
                0228                      globfldtmp2(bi,ip,j,bj,jp) = cbuff(cbuffindex)
                0229 cph)
7109a141b2 Patr*0230 #else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
dff4940422 Patr*0231                      cbuff(cbuffindex) = globfld3d(ii,bi,ip,jj,bj,jp,kk)
7109a141b2 Patr*0232 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
                0233                   endif
                0234                  enddo
                0235                 enddo
                0236                enddo
                0237               enddo
                0238              enddo
                0239 c           --> check cbuffindex.
                0240              if ( cbuffindex .gt. 0) then
                0241                 write(cunit) cbuffindex
                0242                 write(cunit) k
                0243                 write(cunit) (cbuff(ii), ii=1,cbuffindex)
                0244              endif
                0245 c
f5224d0b03 Patr*0246             if ( doPackDiag ) then
                0247                write(cunit2,rec=irectrue) globfldtmp2
                0248                write(cunit3,rec=irectrue) globfldtmp3
                0249             endif
                0250 c
7109a141b2 Patr*0251 c     -- end of k loop --
                0252             enddo
                0253 c     -- end of i loop --
                0254          enddo
                0255 c     -- end of irec loop --
                0256       enddo
                0257 
5cf4364659 Mart*0258       do irec = nrec_nl*nx+1, ncvarrecs(ivar)
5b7a03205a Mart*0259 c     Need to solve for what iobcs would have been.
7109a141b2 Patr*0260          iobcs= mod(irec-1,nobcs)+1
5d5c0b0d52 Patr*0261 
                0262          call MDSREADFIELD_YZ_GL( fname, ctrlprec, 'RL',
de57a2ec4b Mart*0263      &        Nr, globfldyz, irec, mythid)
5d5c0b0d52 Patr*0264 
5cf4364659 Mart*0265          write(cunit) ncvarindex(ivar)
5d5c0b0d52 Patr*0266          write(cunit) 1
                0267          write(cunit) 1
de57a2ec4b Mart*0268          do k = 1,Nr
                0269             irectrue = (irec-1)*nobcs*Nr + (iobcs-1)*Nr + k
5d5c0b0d52 Patr*0270             cbuffindex = 0
                0271             do jp = 1,nPy
                0272              do bj = jtlo,jthi
                0273               do ip = 1,nPx
                0274                do bi = itlo,ithi
                0275                 do j = jmin,jmax
7109a141b2 Patr*0276                  if (globmskyz(bi,ip,j,bj,jp,k,iobcs)  .ne. 0. ) then
5d5c0b0d52 Patr*0277                      cbuffindex = cbuffindex + 1
f5224d0b03 Patr*0278 cph(
                0279                      globfldtmp3(bi,ip,j,bj,jp) =
                0280      &                    globfldyz(bi,ip,j,bj,jp,k)
                0281 cph)
5d5c0b0d52 Patr*0282 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                0283                      if (lxxadxx) then
951926fb9b Jean*0284                         cbuff(cbuffindex) =
5d5c0b0d52 Patr*0285      &                       globfldyz(bi,ip,j,bj,jp,k) *
7109a141b2 Patr*0286 #ifdef CTRL_PACK_PRECISE
                0287      &                       sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
                0288 #else
5d5c0b0d52 Patr*0289      &                       sqrt(weightfld(k,iobcs))
7109a141b2 Patr*0290 #endif
5d5c0b0d52 Patr*0291                      else
951926fb9b Jean*0292                         cbuff(cbuffindex) =
5d5c0b0d52 Patr*0293      &                       globfldyz(bi,ip,j,bj,jp,k) /
7109a141b2 Patr*0294 #ifdef CTRL_PACK_PRECISE
                0295      &                       sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
                0296 #else
5d5c0b0d52 Patr*0297      &                       sqrt(weightfld(k,iobcs))
7109a141b2 Patr*0298 #endif
5d5c0b0d52 Patr*0299                      endif
f5224d0b03 Patr*0300 cph(
                0301                      globfldtmp2(bi,ip,j,bj,jp) = cbuff(cbuffindex)
                0302 cph)
7109a141b2 Patr*0303 #else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
5d5c0b0d52 Patr*0304                      cbuff(cbuffindex) = globfldyz(bi,ip,j,bj,jp,k)
7109a141b2 Patr*0305 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
5d5c0b0d52 Patr*0306                  endif
                0307                 enddo
                0308                enddo
                0309               enddo
                0310              enddo
                0311             enddo
                0312 c           --> check cbuffindex.
                0313             if ( cbuffindex .gt. 0) then
                0314                write(cunit) cbuffindex
                0315                write(cunit) k
                0316                write(cunit) (cbuff(ii), ii=1,cbuffindex)
                0317             endif
                0318 c
f5224d0b03 Patr*0319             if ( doPackDiag ) then
                0320                write(cunit2,rec=irectrue) globfldtmp2
                0321                write(cunit3,rec=irectrue) globfldtmp3
                0322             endif
                0323 c
7109a141b2 Patr*0324 c     -- end of k loop --
                0325          enddo
5d5c0b0d52 Patr*0326 c     -- end of irec loop --
                0327       enddo
                0328 
                0329       _END_MASTER( mythid )
                0330 
dac57cef35 Patr*0331 #endif
                0332 
5d5c0b0d52 Patr*0333       return
                0334       end