Back to home page

MITgcm

 
 

    


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