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_xz(
5cf4364659 Mart*0004      &     cunit, ivar, fname, masktype, weighttype,
5d5c0b0d52 Patr*0005      &     weightfld, nwetglobal, mythid)
                0006 
                0007 c     ==================================================================
                0008 c     SUBROUTINE ctrl_set_unpack_xz
                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
abacc7d9db 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      ( sNx*nSx*nPx*nSy*nPy )
                0066       real*4 globfldtmp2( sNx,nSx,nPx,nSy,nPy )
                0067       real*4 globfldtmp3( sNx,nSx,nPx,nSy,nPy )
                0068       _RL     globfldxz( sNx,nSx,nPx,nSy,nPy,Nr )
                0069       _RL     globfld3d( sNx,nSx,nPx,sNy,nSy,nPy,Nr )
                0070       _RL     globmskxz( sNx,nSx,nPx,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   weightfldxz( sNx,nSx,nPx,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 ip = 1,nPx
                0112                   do bi = itlo,ithi
                0113                      do i = imin,imax
318307105a Patr*0114                         globfldxz  (i,bi,ip,bj,jp,k) = 0. _d 0
                0115                         globfldtmp2(i,bi,ip,bj,jp)   = 0.
                0116                         globfldtmp3(i,bi,ip,bj,jp)   = 0.
7109a141b2 Patr*0117                         do iobcs=1,nobcs
                0118                            globmskxz(i,bi,ip,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 = sNx*nSx*nPx*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
951926fb9b Jean*0173          call MDSREADFIELD_XZ_GL(
5d5c0b0d52 Patr*0174      &        masktype, ctrlprec, 'RL',
7109a141b2 Patr*0175      &        Nr, globmskxz(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_XZ_GL(
                0180      &        weightname, ctrlprec, 'RL',
                0181      &        Nr, weightfldxz(1,1,1,1,1,1,iobcs), iobcs, mythid)
                0182 #endif /* CTRL_UNPACK_PRECISE */
                0183       enddo
                0184 
b94917e7bd Mart*0185       if ( useSingleCPUio ) then
                0186 C     MDSWRITEFIELD_XZ_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)/Ny)
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 j=1,sNy
                0195             iobcs= mod((irec-1)*sNy+j-1,nobcs)+1
951926fb9b Jean*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:xz: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_xz'
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
7109a141b2 Patr*0210                cbuffindex = nwetglobal(k,iobcs)
                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_xz'
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_xz'
7109a141b2 Patr*0223                   endif
                0224                   read(cunit) (cbuff(ii), ii=1,cbuffindex)
                0225                endif
318307105a Patr*0226 
7109a141b2 Patr*0227                cbuffindex = 0
de57a2ec4b Mart*0228                jj=mod((j-1)*Nr+k-1,sNy)+1
                0229                kk=int((j-1)*Nr+k-1)/sNy+1
7109a141b2 Patr*0230                do jp = 1,nPy
                0231                 do bj = jtlo,jthi
                0232                  do ip = 1,nPx
                0233                   do bi = itlo,ithi
                0234                    do i = imin,imax
                0235                     if ( globmskxz(i,bi,ip,bj,jp,k,iobcs) .ne. 0. ) then
                0236                        cbuffindex = cbuffindex + 1
951926fb9b Jean*0237                        globfld3d(i,bi,ip,jj,bj,jp,kk) =
abacc7d9db Patr*0238      &                      cbuff(cbuffindex)
318307105a Patr*0239 cph(
7bfe6112e8 Jean*0240                        globfldtmp2(i,bi,ip,bj,jp) =
318307105a Patr*0241      &                      cbuff(cbuffindex)
                0242 cph)
7109a141b2 Patr*0243 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
951926fb9b Jean*0244                        globfld3d(i,bi,ip,jj,bj,jp,kk) =
abacc7d9db Patr*0245      &                      globfld3d(i,bi,ip,jj,bj,jp,kk)/
7109a141b2 Patr*0246 # ifdef CTRL_UNPACK_PRECISE
                0247      &                      sqrt(weightfldxz(i,bi,ip,bj,jp,k,iobcs))
                0248 # else
                0249      &                      sqrt(weightfld(k,iobcs))
                0250 # endif
                0251 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
                0252                     else
abacc7d9db Patr*0253                        globfld3d(i,bi,ip,jj,bj,jp,kk) = 0. _d 0
7109a141b2 Patr*0254                     endif
318307105a Patr*0255 cph(
                0256                     globfldtmp3(i,bi,ip,bj,jp) =
                0257      &                   globfld3d(i,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 j loop --
                0273          enddo
951926fb9b Jean*0274 
7109a141b2 Patr*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*ny+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:xz: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_xz'
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
5d5c0b0d52 Patr*0299             cbuffindex = nwetglobal(k,iobcs)
                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_xz'
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_xz'
5d5c0b0d52 Patr*0312                endif
                0313                read(cunit) (cbuff(ii), ii=1,cbuffindex)
                0314             endif
318307105a Patr*0315 
5d5c0b0d52 Patr*0316             cbuffindex = 0
                0317             do jp = 1,nPy
                0318              do bj = jtlo,jthi
                0319               do ip = 1,nPx
                0320                do bi = itlo,ithi
                0321                 do i = imin,imax
7109a141b2 Patr*0322                   if ( globmskxz(i,bi,ip,bj,jp,k,iobcs) .ne. 0. ) then
5d5c0b0d52 Patr*0323                      cbuffindex = cbuffindex + 1
                0324                      globfldxz(i,bi,ip,bj,jp,k) = cbuff(cbuffindex)
318307105a Patr*0325 cph(
                0326                      globfldtmp2(i,bi,ip,bj,jp) = cbuff(cbuffindex)
                0327 cph)
5d5c0b0d52 Patr*0328 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
951926fb9b Jean*0329                      globfldxz(i,bi,ip,bj,jp,k) =
5d5c0b0d52 Patr*0330      &                    globfldxz(i,bi,ip,bj,jp,k)/
7109a141b2 Patr*0331 # ifdef CTRL_UNPACK_PRECISE
                0332      &                    sqrt(weightfldxz(i,bi,ip,bj,jp,k,iobcs))
                0333 # else
5d5c0b0d52 Patr*0334      &                    sqrt(weightfld(k,iobcs))
7109a141b2 Patr*0335 # endif
                0336 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
5d5c0b0d52 Patr*0337                   else
                0338                      globfldxz(i,bi,ip,bj,jp,k) = 0. _d 0
                0339                   endif
318307105a Patr*0340 cph(
                0341                   globfldtmp3(i,bi,ip,bj,jp) =
                0342      &                 globfldxz(i,bi,ip,bj,jp,k)
                0343 cph)
5d5c0b0d52 Patr*0344                 enddo
                0345                enddo
                0346               enddo
                0347              enddo
                0348             enddo
                0349 c
318307105a Patr*0350             if ( doPackDiag ) then
                0351                write(cunit2,rec=irectrue) globfldtmp2
                0352                write(cunit3,rec=irectrue) globfldtmp3
                0353             endif
                0354 c
7109a141b2 Patr*0355 c     -- end of k loop --
5d5c0b0d52 Patr*0356          enddo
951926fb9b Jean*0357 
5d5c0b0d52 Patr*0358          call MDSWRITEFIELD_XZ_GL( fname, ctrlprec, 'RL',
                0359      &                             Nr, globfldxz, irec,
                0360      &                             optimcycle, mythid)
                0361 
                0362 c     -- end of irec loop --
                0363       enddo
                0364 
                0365       _END_MASTER( mythid )
                0366 
dac57cef35 Patr*0367 #endif
                0368 
5d5c0b0d52 Patr*0369       return
                0370       end