pkg/profiles/profiles_interp.F

24462d2fa8 Patr*0001 #include "PROFILES_OPTIONS.h"
c9bf163375 Ivan*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
5a362746b1 Gael*0005 #ifdef ALLOW_ECCO
                0006 # include "ECCO_OPTIONS.h"
                0007 #endif
c9bf163375 Ivan*0008 #ifdef ALLOW_BLING
                0009 # include "BLING_OPTIONS.h"
                0010 #endif
6a770e0a24 Patr*0011
c9bf163375 Ivan*0012 CBOP
                0013 C !ROUTINE: PROFILES_INTERP
6a770e0a24 Patr*0014
c9bf163375 Ivan*0015 C !INTERFACE: ==========================================================
                0016       SUBROUTINE PROFILES_INTERP(
6a770e0a24 Patr*0017      O traj_cur_out,
c9bf163375 Ivan*0018      I i_cur, j_cur,
39ce977435 Gael*0019      I weights_cur,
cf16ba6028 Gael*0020      I var_cur,
                0021      I itr_cur,
6a770e0a24 Patr*0022      I file_cur,
13d362b8c1 Ou W*0023      I myTime,
c9bf163375 Ivan*0024      I bi, bj,
                0025      I myThid )
6a770e0a24 Patr*0026
c9bf163375 Ivan*0027 C !DESCRIPTION:
                0028 C 3D interpolation of model counterparts for netcdf profiles data
6a770e0a24 Patr*0029
c9bf163375 Ivan*0030 C !USES: ===============================================================
                0031       IMPLICIT NONE
                0032 C     == Global variables ===
6a770e0a24 Patr*0033 #include "EEPARAMS.h"
                0034 #include "SIZE.h"
c9bf163375 Ivan*0035 #include "PARAMS.h"
6a770e0a24 Patr*0036 #include "GRID.h"
                0037 #include "DYNVARS.h"
5a362746b1 Gael*0038 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0039 # include "ECCO_SIZE.h"
                0040 # include "ECCO.h"
24462d2fa8 Patr*0041 #endif
d28c90138c Patr*0042 #ifdef ALLOW_PROFILES
6328b73337 Gael*0043 # include "PROFILES_SIZE.h"
6e4c90fea3 Patr*0044 # include "profiles.h"
                0045 #endif
7cad2d9e22 Gael*0046 #ifdef ALLOW_PTRACERS
c9bf163375 Ivan*0047 # include "PTRACERS_SIZE.h"
                0048 # include "PTRACERS_FIELDS.h"
                0049 #endif
                0050 #ifdef ALLOW_BLING
                0051 # include "BLING_VARS.h"
7cad2d9e22 Gael*0052 #endif
c9bf163375 Ivan*0053
                0054 C !INPUT PARAMETERS: ===================================================
                0055 C i/j_cur     :: i/j - index arrays of model grid points
                0056 C weights_cur :: weights array for profile data
                0057 C var_cur     :: current variable name
                0058 C itr_cur     :: index of current variable
                0059 C file_cur    :: current filenumber
                0060 C myTime      :: Current time in simulation
                0061 C bi,bj       :: Tile indices
                0062 C myThid      :: my Thread Id number
                0063 #ifdef ALLOW_PROFILES
                0064       INTEGER i_cur(NUM_INTERP_POINTS)
                0065       INTEGER j_cur(NUM_INTERP_POINTS)
                0066       _RL weights_cur(NUM_INTERP_POINTS)
                0067 #else
                0068       INTEGER i_cur, j_cur
                0069       _RL weights_cur
                0070 #endif
                0071       CHARACTER*(8) var_cur
                0072       INTEGER itr_cur
                0073       INTEGER file_cur
13d362b8c1 Ou W*0074       _RL  myTime
c9bf163375 Ivan*0075       INTEGER bi, bj
                0076       INTEGER myThid
                0077
                0078 C !OUTPUT PARAMETERS: ==================================================
                0079 C traj_cur_out :: profile of interpolated model data
                0080 #ifdef ALLOW_PROFILES
                0081       _RL traj_cur_out(NLEVELMAX)
5a362746b1 Gael*0082 #else
c9bf163375 Ivan*0083       _RL traj_cur_out
5a362746b1 Gael*0084 #endif
6e4c90fea3 Patr*0085
c9bf163375 Ivan*0086 #ifdef ALLOW_PROFILES
                0087 C !FUNCTIONS: ==========================================================
                0088       INTEGER ILNBLNK
                0089       EXTERNAL ILNBLNK
                0090
                0091 C !LOCAL VARIABLES: ====================================================
39ce977435 Gael*0092       _RL tab_coeffs1(NUM_INTERP_POINTS)
                0093       _RL tab_coeffs3(NUM_INTERP_POINTS)
c9bf163375 Ivan*0094       _RL ponderations(NUM_INTERP_POINTS), pondsSUM
                0095       INTEGER q, k, kk, kcur, ierr
                0096       _RL traj_cur(Nr), mask_cur(Nr)
39ce977435 Gael*0097       _RL tmp_coeff
c9bf163375 Ivan*0098 CEOP
3c8dcfdea9 Gael*0099
c9bf163375 Ivan*0100       ierr = 0
                0101 C-- horizontal interpolation:
                0102       DO k=1,Nr
                0103        pondsSUM=0. _d 0
                0104        DO q=1,NUM_INTERP_POINTS
                0105         IF (var_cur.EQ.'theta') THEN
                0106          tab_coeffs1(q)=theta(i_cur(q),j_cur(q),k,bi,bj)
                0107         ELSEIF (var_cur.EQ.'salt') THEN
                0108          tab_coeffs1(q)=salt(i_cur(q),j_cur(q),k,bi,bj)
                0109         ELSEIF (var_cur.EQ.'pTracer') THEN
7cad2d9e22 Gael*0110 #ifdef ALLOW_PTRACERS
c9bf163375 Ivan*0111          tab_coeffs1(q)=pTracer(i_cur(q),j_cur(q),k,bi,bj,itr_cur)
7cad2d9e22 Gael*0112 #else
c9bf163375 Ivan*0113          tab_coeffs1(q)=0. _d 0
                0114 #endif
                0115 #ifdef ALLOW_BLING
                0116         ELSEIF (var_cur.EQ.'PCO') THEN
                0117                tab_coeffs1(q)=pCO2(i_cur(q),j_cur(q),bi,bj)
                0118         ELSEIF (var_cur.EQ.'PH') THEN
                0119                tab_coeffs1(q)=pH(i_cur(q),j_cur(q),k,bi,bj)
                0120         ELSEIF (var_cur.EQ.'CHL') THEN
                0121                tab_coeffs1(q)=CHL(i_cur(q),j_cur(q),k,bi,bj)
                0122         ELSEIF (var_cur.EQ.'POC') THEN
                0123                tab_coeffs1(q)=POC(i_cur(q),j_cur(q),k,bi,bj)
7cad2d9e22 Gael*0124 #endif
5a362746b1 Gael*0125 #ifdef ALLOW_ECCO
c9bf163375 Ivan*0126         ELSEIF (var_cur.EQ.'eta') THEN
                0127          tab_coeffs1(q)=m_eta(i_cur(q),j_cur(q),bi,bj)
da0e0765ef An T*0128 #endif
c9bf163375 Ivan*0129         ELSEIF (var_cur.EQ.'UE') THEN
                0130          tab_coeffs1(q)=m_UE(i_cur(q),j_cur(q),k,bi,bj)
                0131         ELSEIF (var_cur.EQ.'VN') THEN
                0132          tab_coeffs1(q)=m_VN(i_cur(q),j_cur(q),k,bi,bj)
                0133         ELSE
                0134          tab_coeffs1(q)=0. _d 0
                0135         ENDIF
39ce977435 Gael*0136         tab_coeffs3(q)=maskC(i_cur(q),j_cur(q),k,bi,bj)
6a770e0a24 Patr*0137
c9bf163375 Ivan*0138         ponderations(q) = tab_coeffs3(q)*weights_cur(q)
                0139         pondsSUM = pondsSUM+ponderations(q)
                0140        ENDDO
f527c11034 Gael*0141
c9bf163375 Ivan*0142        IF (pondsSUM.GT.zeroRL) THEN
                0143         traj_cur(k)=0. _d 0
                0144         mask_cur(k)=1. _d 0
                0145         pondsSUM = 1. _d 0/pondsSUM
                0146         DO q=1,NUM_INTERP_POINTS
                0147          traj_cur(k)=traj_cur(k)
                0148      &        +tab_coeffs1(q)*ponderations(q)*pondsSUM
                0149         ENDDO
                0150        ELSE
                0151         traj_cur(k)=0. _d 0
                0152         mask_cur(k)=0. _d 0
                0153        ENDIF
                0154       ENDDO
6a770e0a24 Patr*0155
c9bf163375 Ivan*0156 #ifdef ALLOW_AUTODIFF_TAMC
                0157 C     This is necessary, because TAF does not recognise prof_mask1D_cur
                0158 C     as active and does not include it in the TL-version of this
                0159 C     routine leading to prof_mask1D_cur=0 and also to
                0160 C     objf_profiles_tl=0.
                0161 C$TAF INCOMPLETE prof_mask1D_cur
                0162 #endif
                0163 C--   vertical interpolation:
                0164       DO kk=1,NLEVELMAX
                0165        traj_cur_out(kk)=0. _d 0
                0166        prof_mask1D_cur(kk,bi,bj)=0. _d 0
                0167       ENDDO
                0168       DO kk=1,ProfDepthNo(file_cur,bi,bj)
                0169 C--   case 1: above first grid center=> first grid center value
                0170        IF (prof_depth(file_cur,kk,bi,bj).LT.-rC(1)) THEN
                0171         traj_cur_out(kk)=traj_cur(1)
                0172         prof_mask1D_cur(kk,bi,bj)=mask_cur(1)
                0173 C--   case 2: just below last grid center=> last cell value
                0174        ELSEIF (prof_depth(file_cur,kk,bi,bj).GE.-rC(Nr)) THEN
                0175         IF ( prof_depth(file_cur,kk,bi,bj) .LT.
                0176      &       (-rC(Nr)+drC(Nr)*halfRS) ) THEN
                0177          traj_cur_out(kk)=traj_cur(Nr)
                0178          prof_mask1D_cur(kk,bi,bj)=mask_cur(Nr)
                0179         ENDIF
                0180 C--   case 3: between two grid centers
                0181        ELSE
                0182         kcur=0
                0183         DO k=1,Nr-1
                0184          IF ( (prof_depth(file_cur,kk,bi,bj).GE.-rC(k)).AND.
                0185      &        (prof_depth(file_cur,kk,bi,bj).LT.-rC(k+1))) THEN
                0186           kcur=k
                0187          ENDIF
                0188         ENDDO
                0189         IF (kcur.EQ.0) THEN
                0190          WRITE(errorMessageUnit,'(A)')
b00d6c1700 Gael*0191      &        'ERROR in PROFILES_INTERP: unexpected case 1'
c9bf163375 Ivan*0192          ierr = ierr + 1
                0193         ENDIF
                0194         IF (mask_cur(kcur+1).EQ.oneRL) THEN
                0195 C--   subcase 1: 2 wet points=>linear interpolation
                0196          tmp_coeff=(prof_depth(file_cur,kk,bi,bj)+rC(kcur))/
                0197      &        (-rC(kcur+1)+rC(kcur))
                0198          traj_cur_out(kk)=(1-tmp_coeff)*traj_cur(kcur)
                0199      &        +tmp_coeff*traj_cur(kcur+1)
                0200          prof_mask1D_cur(kk,bi,bj)=1. _d 0
                0201          IF (mask_cur(kcur).EQ.zeroRL) THEN
                0202           WRITE(errorMessageUnit,'(A)')
                0203      &         'ERROR in PROFILES_INTERP: unexpected case 2'
                0204           ierr = ierr + 1
                0205          ENDIF
                0206         ELSEIF (prof_depth(file_cur,kk,bi,bj).LT.-rF(kcur+1)) THEN
                0207 C--   subcase 2: only 1 wet point just above=>upper cell value
                0208          traj_cur_out(kk)=traj_cur(kcur)
                0209          prof_mask1D_cur(kk,bi,bj)=mask_cur(kcur)
                0210         ENDIF
                0211        ENDIF
                0212       ENDDO
                0213
                0214       IF ( ierr .GT. 0 ) THEN
                0215 C--    cannot call ALL_PROC_DIE from here: see issue #439, 2021-03-05 comment
                0216 c      CALL ALL_PROC_DIE( myThid )
                0217        STOP 'ABNORMAL END: S/R PROFILES_INTERP'
                0218       ENDIF
                0219 #endif /* ALLOW_PROFILES */
6a770e0a24 Patr*0220
af9e0d08c7 Jean*0221       RETURN
                0222       END