pkg/layers/layers_init_fixed.F

108a00eab9 Ryan*0001 #include "LAYERS_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004
                0005       SUBROUTINE LAYERS_INIT_FIXED( myThid )
                0006
                0007 C ===================================================================
                0008 C     Initialize LAYERS variables that are kept fixed during the run.
                0009 C ===================================================================
                0010
                0011       IMPLICIT NONE
                0012 #include "EEPARAMS.h"
                0013 #include "SIZE.h"
                0014 #include "PARAMS.h"
                0015 #include "GRID.h"
                0016 #include "LAYERS_SIZE.h"
                0017 #include "LAYERS.h"
                0018
                0019 C  INPUT/OUTPUT PARAMETERS:
                0020 C     myThid ::  my Thread Id number
                0021       INTEGER myThid
                0022
                0023 C  LOCAL VARIABLES:
                0024 C     k         :: loop index
                0025 C     kk,kkinit :: fine grid loop indices
                0026 C     Zf        :: depth at cell boundaries
                0027 C     Zf        :: depth at cell centers
                0028 C     ZZf       :: depth at cell boundaries (fine grid)
                0029 C     ZZc       :: depth at cell centers (fine grid)
14cf1d4767 Jean*0030 C     msgBuf    :: Informational/error message buffer
406891c1a3 Gael*0031       INTEGER k,kk,kkinit,iLa
108a00eab9 Ryan*0032       _RL     Zf(Nr+1)
                0033       _RL     Zc(Nr)
                0034       _RL     ZZf(FineGridMax+1)
                0035       _RL     ZZc(FineGridMax)
                0036
f0f170d54b Mart*0037       CHARACTER*11   tmpName
108a00eab9 Ryan*0038       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0039
406891c1a3 Gael*0040 C     Functions ::
                0041       INTEGER     ILNBLNK
                0042       EXTERNAL    ILNBLNK
f0f170d54b Mart*0043
                0044 #if defined ALLOW_MNC && defined LAYERS_MNC
108a00eab9 Ryan*0045       IF (layers_MNC) THEN
                0046         CALL LAYERS_MNC_INIT( myThid )
                0047       ENDIF
f0f170d54b Mart*0048 #endif
108a00eab9 Ryan*0049
                0050 C  Set up the vertical grid
                0051
                0052 C     for now, just use up the entire available array for ZZ
                0053       NZZ = FineGridMax
                0054
                0055 C     Z and ZZ are INCREASING DOWNWARD!!!
                0056 C     Maybe this is dumb but it will work as long as we are consistent
                0057
                0058 C     Each dF cell is split into FineGridFact fine cells
                0059 C     Calculate dZZf on the fine grid
                0060       kkinit = 1
                0061       DO k=1,Nr
                0062         DO kk=kkinit,kkinit+FineGridFact-1
                0063           dZZf(kk) = dRf(k) / FineGridFact
                0064         ENDDO
                0065         kkinit = kkinit + FineGridFact
                0066       ENDDO
                0067
                0068 C     find the depths
                0069       Zf(1) = 0. _d 0
                0070       Zc(1) = drC(1)
                0071       DO k=2,Nr
                0072         Zf(k) = Zf(k-1) + drF(k-1)
                0073         Zc(k) = Zc(k-1) + drC(k)
                0074       ENDDO
                0075       Zf(Nr+1) = Zf(Nr) + drF(Nr)
                0076
                0077 C     create ZZ
                0078       ZZf(1) = 0. _d 0
f61c9be03a Ryan*0079       ZZc(1) = 0.5 _d 0 * dZZf(1)
108a00eab9 Ryan*0080
f61c9be03a Ryan*0081       DO kk=2,NZZ+1
108a00eab9 Ryan*0082             ZZf(kk) = ZZf(kk-1) + dZZf(kk-1)
f61c9be03a Ryan*0083             ZZc(kk-1) = 0.5 _d 0 * (ZZf(kk) + ZZf(kk-1))
108a00eab9 Ryan*0084       ENDDO
                0085
                0086 C     create the interpolating mapping arrays
                0087       k = 1
                0088       DO kk=1,NZZ
                0089 C       see if ZZc point is less than the top Zc point
                0090         IF ( ZZc(kk) .LT. Zc(1) ) THEN
                0091           MapIndex(kk) = 1
                0092           MapFact(kk) = 1.0 _d 0
                0093 C       see if ZZc point is greater than the bottom Zc point
f0f170d54b Mart*0094         ELSEIF ( (ZZc(kk) .GE. Zc(Nr)) .OR. (k .EQ. Nr) ) THEN
108a00eab9 Ryan*0095           MapIndex(kk) = Nr - 1
                0096           MapFact(kk) = 0.0 _d 0
                0097 C       Otherwise we are somewhere in between and need to do interpolation)
f0f170d54b Mart*0098         ELSEIF ( (ZZc(kk) .GE. Zc(k))
108a00eab9 Ryan*0099      &   .AND. (ZZc(kk) .LT. Zc(Nr)) ) THEN
                0100 C         Find the proper k value
                0101           DO WHILE (ZZc(kk) .GE. Zc(k+1))
                0102             k = k + 1
                0103           ENDDO
                0104 C         If the loop stopped, that means Zc(k) <= ZZc(kk) < ZZc(k+1)
                0105           MapIndex(kk) = k
                0106           MapFact(kk) = 1.0 - (( ZZc(kk) - Zc(k) ) / drC(k+1))
                0107         ELSE
                0108 C       This means there was a problem
                0109           WRITE(msgBuf,'(A,I4,A,I4,A,1E14.6,A,2E14.6)')
                0110      &     'S/R LAYERS_INIT_FIXED: kk=', kk, ', k=', k,
                0111      &     ', ZZc(kk)=', ZZc(kk),' , Zc(k)=',Zc(k)
                0112           CALL PRINT_ERROR( msgBuf, myThid )
                0113           STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'
14cf1d4767 Jean*0114         ENDIF
108a00eab9 Ryan*0115
                0116 C       See which grid box the point lies in
e643b51506 Jean*0117         IF ( ZZc(kk).LT.Zf(MapIndex(kk)+1) ) THEN
108a00eab9 Ryan*0118           CellIndex(kk) = MapIndex(kk)
                0119         ELSE
                0120           CellIndex(kk) = MapIndex(kk)+1
14cf1d4767 Jean*0121         ENDIF
108a00eab9 Ryan*0122       ENDDO
                0123
8830b8f970 Jean*0124       IF ( debugLevel .GE. debLevB ) THEN
108a00eab9 Ryan*0125         CALL PRINT_MESSAGE( 'LAYERS_INIT_FIXED Debugging:',
8830b8f970 Jean*0126      &             standardMessageUnit, SQUEEZE_RIGHT, myThid )
108a00eab9 Ryan*0127         DO kk=1,NZZ
                0128           WRITE(msgBuf,'(A,1F6.1,A,I3,A,I3,A,I3,A,1F6.4,A,I3,A,I3)')
                0129      &     '// ZZc=', ZZc(kk),
                0130      &     ', MapIndex(',kk,')=',MapIndex(kk),
                0131      &     ', MapFact(',kk,')=',MapFact(kk),
                0132      &     ', CellIndex(',kk,')=',CellIndex(kk)
                0133           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
8830b8f970 Jean*0134      &                        SQUEEZE_RIGHT, myThid )
108a00eab9 Ryan*0135         ENDDO
14cf1d4767 Jean*0136       ENDIF
108a00eab9 Ryan*0137
f61c9be03a Ryan*0138 C     Output the layer coordinates
406891c1a3 Gael*0139       DO iLa=1,layers_maxNum
                0140       IF ( layers_num(iLa).NE.0 ) THEN
                0141       WRITE(tmpName,'(A7,I1,A3)') 'layers',iLa,layers_name(iLa)
f8fe443de1 Davi*0142       CALL WRITE_GLVEC_RL( tmpName, ' ',layers_bounds(1,iLa),1+Nlayers,
f61c9be03a Ryan*0143      & -1, myThid )
406891c1a3 Gael*0144       ENDIF
                0145       ENDDO
f0f170d54b Mart*0146
                0147 C--   Set up layers "w-grid" for transformation calculation
50d8304171 Ryan*0148 #ifdef LAYERS_THERMODYNAMICS
                0149       DO iLa=1,layers_maxNum
                0150        IF ( layers_num(iLa).NE.0 ) THEN
                0151         DO k=1,Nlayers
                0152           layers_bounds_w(k,iLa) = 0.5 _d 0 * (
                0153      &           layers_bounds(k+1,iLa) +
                0154      &           layers_bounds(k,iLa) )
                0155         ENDDO
                0156         DO k=1,Nlayers-1
                0157           layers_recip_delta(k,iLa) = 1.0 _d 0 / (
                0158      &           layers_bounds_w(k+1,iLa) -
a3d5c99c9c Ryan*0159      &           layers_bounds_w(k,iLa) )
50d8304171 Ryan*0160         ENDDO
                0161        ENDIF
                0162       ENDDO
                0163 #endif /* LAYERS_THERMODYNAMICS */
f61c9be03a Ryan*0164
6e34151e15 Gael*0165 #ifdef ALLOW_DIAGNOSTICS
                0166       IF ( useDiagnostics ) THEN
                0167         CALL LAYERS_DIAGNOSTICS_INIT( myThid )
                0168       ENDIF
                0169 #endif
                0170
108a00eab9 Ryan*0171       RETURN
                0172       END