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