pkg/ctrl/ctrl_set_pack_xy.F

7bfe6112e8 Jean*0001 #include "CTRL_OPTIONS.h"
5d5c0b0d52 Patr*0002
                0003       subroutine ctrl_set_pack_xy(
5cf4364659 Mart*0004      &     cunit, ivar, precondScalar,
3145d51f22 Patr*0005      &     fname, masktype, weighttype,
4d72283393 Mart*0006      &     lxxadxx, myThid )
5d5c0b0d52 Patr*0007
                0008 c     ==================================================================
                0009 c     SUBROUTINE ctrl_set_pack_xy
                0010 c     ==================================================================
                0011 c
                0012 c     o Compress the control vector such that only ocean points are
                0013 c       written to file.
                0014 c
                0015 c     ==================================================================
                0016
                0017       implicit none
                0018
                0019 c     == global variables ==
                0020
                0021 #include "EEPARAMS.h"
                0022 #include "SIZE.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
                0025
5cf4364659 Mart*0026 #include "CTRL_SIZE.h"
4d72283393 Mart*0027 #include "CTRL.h"
65754df434 Mart*0028 #include "OPTIMCYCLE.h"
5d5c0b0d52 Patr*0029
                0030 c     == routine arguments ==
                0031
                0032       integer cunit
5cf4364659 Mart*0033       integer ivar
3145d51f22 Patr*0034       _RL precondScalar
de57a2ec4b Mart*0035       character*(MAX_LEN_FNAM) fname
45913d6a59 Patr*0036       character*(  9) masktype
de57a2ec4b Mart*0037       character*(MAX_LEN_FNAM) weighttype
5d5c0b0d52 Patr*0038       logical lxxadxx
4d72283393 Mart*0039       integer myThid
5d5c0b0d52 Patr*0040
dac57cef35 Patr*0041 #ifndef EXCLUDE_CTRL_PACK
9f5240b52a Jean*0042 c     == external ==
                0043       integer  ilnblnk
                0044       external ilnblnk
5d5c0b0d52 Patr*0045
9f5240b52a Jean*0046 c     == local variables ==
5d5c0b0d52 Patr*0047       integer bi,bj
                0048       integer i,j,k
9f5240b52a Jean*0049       integer ii, irec
                0050       integer cbuffindex
                0051       real*4 cbuff( sNx*nSx*nPx*sNy*nSy*nPy )
de57a2ec4b Mart*0052       character*(MAX_LEN_FNAM) cfile2, cfile3
9f5240b52a Jean*0053 C========================================================================
                0054 # ifndef ALLOW_PACKUNPACK_METHOD2
                0055       integer ip,jp
                0056       integer nrec_nl
5d5c0b0d52 Patr*0057       integer itlo,ithi
                0058       integer jtlo,jthi
                0059       integer jmin,jmax
                0060       integer imin,imax
9f5240b52a Jean*0061       _RL     globmsk  ( sNx,nSx,nPx,sNy,nSy,nPy,Nr )
                0062       _RL     globfld3d( sNx,nSx,nPx,sNy,nSy,nPy,Nr )
0ba65c94ff Patr*0063       integer reclen, irectrue
                0064       integer cunit2, cunit3
9f5240b52a Jean*0065       real*4 globfldtmp2( sNx,nSx,nPx,sNy,nSy,nPy )
                0066       real*4 globfldtmp3( sNx,nSx,nPx,sNy,nSy,nPy )
                0067 # else /* ALLOW_PACKUNPACK_METHOD2 */
                0068       integer il
                0069       _RL msk3d     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0070       real*8 msk2d_buf    (sNx,sNy,nSx,nSy)
                0071       real*8 msk2d_buf_glo(Nx,Ny)
                0072       real*8 fld2d_buf    (sNx,sNy,nSx,nSy)
                0073       real*8 fld2d_buf_glo(Nx,Ny)
                0074       _RL fld2dDim  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0075       _RL fld2dNodim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0076       _RL dummy
                0077 # endif /* ALLOW_PACKUNPACK_METHOD2 */
5d5c0b0d52 Patr*0078 c     == end of interface ==
                0079
9f5240b52a Jean*0080 # ifndef ALLOW_PACKUNPACK_METHOD2
5d5c0b0d52 Patr*0081       jtlo = 1
9f5240b52a Jean*0082       jthi = nSy
5d5c0b0d52 Patr*0083       itlo = 1
9f5240b52a Jean*0084       ithi = nSx
5d5c0b0d52 Patr*0085       jmin = 1
9f5240b52a Jean*0086       jmax = sNy
5d5c0b0d52 Patr*0087       imin = 1
9f5240b52a Jean*0088       imax = sNx
5d5c0b0d52 Patr*0089
8f0b59c61c Patr*0090       nbuffglobal = nbuffglobal + 1
                0091
5d5c0b0d52 Patr*0092 c     Initialise temporary file
9f5240b52a Jean*0093       do k = 1,Nr
3145d51f22 Patr*0094        do jp = 1,nPy
                0095         do bj = jtlo,jthi
                0096          do j = jmin,jmax
                0097           do ip = 1,nPx
                0098            do bi = itlo,ithi
                0099             do i = imin,imax
                0100              globfld3d  (i,bi,ip,j,bj,jp,k) = 0. _d 0
                0101              globmsk    (i,bi,ip,j,bj,jp,k) = 0. _d 0
                0102              globfldtmp2(i,bi,ip,j,bj,jp)   = 0. _d 0
                0103              globfldtmp3(i,bi,ip,j,bj,jp)   = 0. _d 0
5d5c0b0d52 Patr*0104             enddo
3145d51f22 Patr*0105            enddo
                0106           enddo
5d5c0b0d52 Patr*0107          enddo
3145d51f22 Patr*0108         enddo
                0109        enddo
5d5c0b0d52 Patr*0110       enddo
                0111
                0112 c--   Only the master thread will do I/O.
4d72283393 Mart*0113       _BEGIN_MASTER( myThid )
5d5c0b0d52 Patr*0114
0ba65c94ff Patr*0115       if ( doPackDiag ) then
                0116          if ( lxxadxx ) then
de57a2ec4b Mart*0117             write(cfile2,'(a,I3.3,a,I4.4,a)')
951926fb9b Jean*0118      &           'diag_pack_nonout_ctrl_',
5cf4364659 Mart*0119      &           ivar, '_', optimcycle, '.bin'
de57a2ec4b Mart*0120             write(cfile3,'(a,I3.3,a,I4.4,a)')
951926fb9b Jean*0121      &           'diag_pack_dimout_ctrl_',
5cf4364659 Mart*0122      &           ivar, '_', optimcycle, '.bin'
0ba65c94ff Patr*0123          else
de57a2ec4b Mart*0124             write(cfile2,'(a,I3.3,a,I4.4,a)')
951926fb9b Jean*0125      &           'diag_pack_nonout_grad_',
5cf4364659 Mart*0126      &           ivar, '_', optimcycle, '.bin'
de57a2ec4b Mart*0127             write(cfile3,'(a,I3.3,a,I4.4,a)')
951926fb9b Jean*0128      &           'diag_pack_dimout_grad_',
5cf4364659 Mart*0129      &           ivar, '_', optimcycle, '.bin'
0ba65c94ff Patr*0130          endif
                0131
9f5240b52a Jean*0132          reclen = FLOAT(sNx*nSx*nPx*sNy*nSy*nPy*4)
4d72283393 Mart*0133          call mdsfindunit( cunit2, myThid )
0ba65c94ff Patr*0134          open( cunit2, file=cfile2, status='unknown',
                0135      &        access='direct', recl=reclen )
4d72283393 Mart*0136          call mdsfindunit( cunit3, myThid )
0ba65c94ff Patr*0137          open( cunit3, file=cfile3, status='unknown',
                0138      &        access='direct', recl=reclen )
                0139       endif
                0140
951926fb9b Jean*0141       call MDSREADFIELD_3D_GL(
5d5c0b0d52 Patr*0142      &     masktype, ctrlprec, 'RL',
4d72283393 Mart*0143      &     Nr, globmsk, 1, myThid)
5d5c0b0d52 Patr*0144
5cf4364659 Mart*0145       nrec_nl=int(ncvarrecs(ivar)/Nr)
7109a141b2 Patr*0146       do irec = 1, nrec_nl
                0147
cf705a6c8e Mart*0148        call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
4d72283393 Mart*0149      &      Nr, globfld3d, irec, myThid)
cf705a6c8e Mart*0150        do k = 1, Nr
                0151         irectrue = (irec-1)*Nr + k
faf44775ba Patr*0152 #ifndef ALLOW_ADMTLM
5cf4364659 Mart*0153         write(cunit) ncvarindex(ivar)
cf705a6c8e Mart*0154         write(cunit) 1
                0155         write(cunit) 1
faf44775ba Patr*0156 #endif
cf705a6c8e Mart*0157         cbuffindex = 0
                0158         do jp = 1,nPy
                0159          do bj = jtlo,jthi
                0160           do j = jmin,jmax
                0161            do ip = 1,nPx
                0162             do bi = itlo,ithi
                0163              do i = imin,imax
                0164               if (globmsk(i,bi,ip,j,bj,jp,1)  .ne. 0. ) then
                0165                cbuffindex = cbuffindex + 1
0ba65c94ff Patr*0166 cph(
cf705a6c8e Mart*0167                globfldtmp3(i,bi,ip,j,bj,jp) =
                0168      &              globfld3d(i,bi,ip,j,bj,jp,k)
0ba65c94ff Patr*0169 cph)
cf705a6c8e Mart*0170                if (lxxadxx) then
                0171                 cbuff(cbuffindex) =
                0172      &               globfld3d(i,bi,ip,j,bj,jp,k)
                0173      &               / PrecondScalar
                0174                else
                0175                 cbuff(cbuffindex) =
                0176      &               globfld3d(i,bi,ip,j,bj,jp,k)
                0177      &               * PrecondScalar
                0178                endif
0ba65c94ff Patr*0179 cph(
cf705a6c8e Mart*0180                globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
0ba65c94ff Patr*0181 cph)
faf44775ba Patr*0182 #ifdef ALLOW_ADMTLM
cf705a6c8e Mart*0183                nveccount = nveccount + 1
                0184                phtmpadmtlm(nveccount) = cbuff(cbuffindex)
faf44775ba Patr*0185 #endif
cf705a6c8e Mart*0186               endif
7109a141b2 Patr*0187              enddo
                0188             enddo
cf705a6c8e Mart*0189            enddo
                0190           enddo
                0191          enddo
                0192         enddo
7109a141b2 Patr*0193 c           --> check cbuffindex.
cf705a6c8e Mart*0194         if ( cbuffindex .gt. 0) then
faf44775ba Patr*0195 #ifndef ALLOW_ADMTLM
cf705a6c8e Mart*0196          write(cunit) cbuffindex
                0197          write(cunit) 1
e324d8bc16 Patr*0198 cph#endif
cf705a6c8e Mart*0199          write(cunit) (cbuff(ii), ii=1,cbuffindex)
e324d8bc16 Patr*0200 #endif
cf705a6c8e Mart*0201         endif
0ba65c94ff Patr*0202 c
cf705a6c8e Mart*0203         if ( doPackDiag ) then
                0204          write(cunit2,rec=irectrue) globfldtmp2
                0205          write(cunit3,rec=irectrue) globfldtmp3
                0206         endif
0ba65c94ff Patr*0207 c
cf705a6c8e Mart*0208        enddo
7109a141b2 Patr*0209 c
                0210 c     -- end of irec loop --
                0211       enddo
                0212
5cf4364659 Mart*0213       do irec = nrec_nl*Nr+1, ncvarrecs(ivar)
5d5c0b0d52 Patr*0214
cf705a6c8e Mart*0215        call MDSREADFIELD_2D_GL( fname, ctrlprec, 'RL',
4d72283393 Mart*0216      &        1, globfld3d(1,1,1,1,1,1,1), irec, myThid)
5d5c0b0d52 Patr*0217
faf44775ba Patr*0218 #ifndef ALLOW_ADMTLM
5cf4364659 Mart*0219        write(cunit) ncvarindex(ivar)
cf705a6c8e Mart*0220        write(cunit) 1
                0221        write(cunit) 1
faf44775ba Patr*0222 #endif
cf705a6c8e Mart*0223        do k = 1, 1
                0224         irectrue = irec
                0225         cbuffindex = 0
                0226         do jp = 1,nPy
                0227          do bj = jtlo,jthi
                0228           do j = jmin,jmax
                0229            do ip = 1,nPx
                0230             do bi = itlo,ithi
                0231              do i = imin,imax
                0232               if (globmsk(i,bi,ip,j,bj,jp,k)  .ne. 0. ) then
                0233                cbuffindex = cbuffindex + 1
0ba65c94ff Patr*0234 cph(
cf705a6c8e Mart*0235                globfldtmp3(i,bi,ip,j,bj,jp) =
                0236      &              globfld3d(i,bi,ip,j,bj,jp,k)
0ba65c94ff Patr*0237 cph)
cf705a6c8e Mart*0238                if (lxxadxx) then
                0239                 cbuff(cbuffindex) =
                0240      &               globfld3d(i,bi,ip,j,bj,jp,k)
                0241      &               / PrecondScalar
                0242                else
                0243                 cbuff(cbuffindex) =
                0244      &               globfld3d(i,bi,ip,j,bj,jp,k)
                0245      &               * PrecondScalar
                0246                endif
0ba65c94ff Patr*0247 cph(
cf705a6c8e Mart*0248                globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
0ba65c94ff Patr*0249 cph)
8f0b59c61c Patr*0250 #ifdef ALLOW_ADMTLM
cf705a6c8e Mart*0251                nveccount = nveccount + 1
                0252                phtmpadmtlm(nveccount) = cbuff(cbuffindex)
8f0b59c61c Patr*0253 #endif
cf705a6c8e Mart*0254               endif
5d5c0b0d52 Patr*0255              enddo
                0256             enddo
cf705a6c8e Mart*0257            enddo
                0258           enddo
                0259          enddo
                0260         enddo
5d5c0b0d52 Patr*0261 c           --> check cbuffindex.
                0262             if ( cbuffindex .gt. 0) then
faf44775ba Patr*0263 #ifndef ALLOW_ADMTLM
5d5c0b0d52 Patr*0264                write(cunit) cbuffindex
                0265                write(cunit) k
e324d8bc16 Patr*0266 cph#endif
5d5c0b0d52 Patr*0267                write(cunit) (cbuff(ii), ii=1,cbuffindex)
e324d8bc16 Patr*0268 #endif
5d5c0b0d52 Patr*0269             endif
0ba65c94ff Patr*0270 c
                0271             if ( doPackDiag ) then
                0272                write(cunit2,rec=irectrue) globfldtmp2
                0273                write(cunit3,rec=irectrue) globfldtmp3
                0274             endif
                0275 c
5d5c0b0d52 Patr*0276          enddo
                0277 c
                0278 c     -- end of irec loop --
                0279       enddo
                0280
0ba65c94ff Patr*0281       if ( doPackDiag ) then
                0282          close ( cunit2 )
                0283          close ( cunit3 )
                0284       endif
951926fb9b Jean*0285
4d72283393 Mart*0286       _END_MASTER( myThid )
5d5c0b0d52 Patr*0287
9f5240b52a Jean*0288 # else /* ALLOW_PACKUNPACK_METHOD2 */
23a37235f2 Gael*0289
                0290 c-- part 1: preliminary reads and definitions
                0291
1c8d09be4c Gael*0292 #ifdef ALLOW_AUTODIFF
cf2ce61250 Jean*0293       call active_read_xyz(masktype, msk3d, 1,
4d72283393 Mart*0294      &    .FALSE., .FALSE., 0 , myThid, dummy)
1c8d09be4c Gael*0295 #else
                0296       CALL READ_REC_XYZ_RL( masktype, msk3d, 1, 1, myThid )
                0297 #endif
23a37235f2 Gael*0298
                0299       if ( doPackDiag ) then
                0300          il = ilnblnk( fname )
                0301          if ( lxxadxx ) then
de57a2ec4b Mart*0302             write(cfile2,'(2a)') fname(1:il),'.pack_ctrl_adim'
                0303             write(cfile3,'(2a)') fname(1:il),'.pack_ctrl_dim'
23a37235f2 Gael*0304          else
de57a2ec4b Mart*0305             write(cfile2,'(2a)') fname(1:il),'.pack_grad_adim'
                0306             write(cfile3,'(2a)') fname(1:il),'.pack_grad_dim'
23a37235f2 Gael*0307          endif
                0308       endif
                0309
                0310 c-- part 2: loop over records
                0311
5cf4364659 Mart*0312       do irec = 1, ncvarrecs(ivar)
23a37235f2 Gael*0313
cf2ce61250 Jean*0314 c-- 2.1:
                0315       call READ_REC_3D_RL( fname, ctrlprec,
4d72283393 Mart*0316      &        1, fld2dDim, irec, 0, myThid)
23a37235f2 Gael*0317
                0318 c-- 2.2: normalize field if needed
                0319       DO bj = myByLo(myThid), myByHi(myThid)
                0320        DO bi = myBxLo(myThid), myBxHi(myThid)
                0321          DO j=1,sNy
                0322           DO i=1,sNx
                0323            if (msk3d(i,j,1,bi,bj).EQ.0. _d 0) then
                0324             fld2dDim(i,j,bi,bj)=0. _d 0
                0325             fld2dNodim(i,j,bi,bj)=0. _d 0
                0326            else
                0327             if (lxxadxx) then
cf705a6c8e Mart*0328              fld2dNodim(i,j,bi,bj) =
                0329      &            fld2dDim(i,j,bi,bj) / PrecondScalar
23a37235f2 Gael*0330             else
cf705a6c8e Mart*0331              fld2dNodim(i,j,bi,bj) =
3145d51f22 Patr*0332      &              fld2dDim(i,j,bi,bj) * PrecondScalar
96c2862ec2 Gael*0333             endif
23a37235f2 Gael*0334            endif
                0335           ENDDO
                0336          ENDDO
                0337        ENDDO
                0338       ENDDO
                0339
                0340 c-- 2.3:
                0341       if ( doPackDiag ) then
                0342       call WRITE_REC_3D_RL( cfile2, ctrlprec,
4d72283393 Mart*0343      &        1, fld2dNodim, irec, 0, myThid)
23a37235f2 Gael*0344       call WRITE_REC_3D_RL( cfile3, ctrlprec,
4d72283393 Mart*0345      &        1, fld2dDim, irec, 0, myThid)
23a37235f2 Gael*0346       endif
                0347
                0348 c-- 2.4: array -> buffer -> global buffer -> global file
                0349
                0350 #ifndef ALLOW_ADMTLM
4d72283393 Mart*0351       _BEGIN_MASTER( myThid )
23a37235f2 Gael*0352       IF ( myProcId .eq. 0 ) THEN
5cf4364659 Mart*0353          write(cunit) ncvarindex(ivar)
23a37235f2 Gael*0354          write(cunit) 1
                0355          write(cunit) 1
                0356       ENDIF
4d72283393 Mart*0357       _END_MASTER( myThid )
23a37235f2 Gael*0358       _BARRIER
                0359 #endif
                0360
                0361       do k = 1, 1
                0362
                0363         CALL MDS_PASS_R8toRL( fld2d_buf, fld2dNodim,
cf2ce61250 Jean*0364      &                 0, 0, 1, k, 1, 0, 0, .FALSE., myThid )
23a37235f2 Gael*0365         CALL BAR2( myThid )
                0366         CALL GATHER_2D_R8( fld2d_buf_glo, fld2d_buf,
                0367      &                       Nx,Ny,.FALSE.,.TRUE.,myThid)
                0368         CALL BAR2( myThid )
                0369
                0370         CALL MDS_PASS_R8toRL( msk2d_buf, msk3d,
cf2ce61250 Jean*0371      &                 0, 0, 1, k, Nr, 0, 0, .FALSE., myThid )
23a37235f2 Gael*0372         CALL BAR2( myThid )
                0373         CALL GATHER_2D_R8( msk2d_buf_glo, msk2d_buf,
                0374      &                       Nx,Ny,.FALSE.,.TRUE.,myThid)
                0375         CALL BAR2( myThid )
                0376
4d72283393 Mart*0377         _BEGIN_MASTER( myThid )
23a37235f2 Gael*0378         cbuffindex = 0
                0379         IF ( myProcId .eq. 0 ) THEN
                0380
                0381         DO j=1,Ny
                0382           DO i=1,Nx
                0383             if (msk2d_buf_glo(i,j) .ne. 0. ) then
                0384                cbuffindex = cbuffindex + 1
                0385                cbuff(cbuffindex) = fld2d_buf_glo(i,j)
                0386 #ifdef ALLOW_ADMTLM
                0387                nveccount = nveccount + 1
                0388                phtmpadmtlm(nveccount) = cbuff(cbuffindex)
dac57cef35 Patr*0389 #endif
23a37235f2 Gael*0390             endif
                0391           ENDDO
                0392         ENDDO
                0393
                0394 #ifndef ALLOW_ADMTLM
                0395         if ( cbuffindex .gt. 0) then
                0396           write(cunit) cbuffindex
                0397           write(cunit) k
                0398           write(cunit) (cbuff(ii), ii=1,cbuffindex)
                0399         endif
                0400 #endif
                0401
                0402         ENDIF
4d72283393 Mart*0403         _END_MASTER( myThid )
23a37235f2 Gael*0404         _BARRIER
                0405
                0406       enddo
                0407       enddo
dac57cef35 Patr*0408
cf2ce61250 Jean*0409 # endif /* ALLOW_PACKUNPACK_METHOD2 */
23a37235f2 Gael*0410 # endif /* EXCLUDE_CTRL_PACK */
cf2ce61250 Jean*0411
5d5c0b0d52 Patr*0412       return
                0413       end