pkg/ctrl/ctrl_set_pack_xyz.F

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