Allow using oasim pkg with radtrans

pkg/radtrans/RADTRANS_PARAMS.h

@@ -108,7 +108,8 @@ C                       if false use angle at actual time
      &    RT_wbTotalWidth,
      &    RT_WtouEins,
      &    RT_rmus, RT_rmuu, RT_rd, RT_ru,
-     &    RT_sfcIrrThresh
+     &    RT_sfcIrrThresh,
+     &    RT_oasimWgt
       _RL RT_refract_water
       _RL RT_rmud_max
       _RL RT_wbRefWLs(nlam)
@@ -118,15 +119,18 @@ C                       if false use angle at actual time
       _RL RT_WtouEins(nlam)
       _RL RT_rmus, RT_rmuu, RT_rd, RT_ru
       _RL RT_sfcIrrThresh
+      _RL RT_oasimWgt(nlam)
 
       COMMON/radtrans_params_i/ RT_kmax
       INTEGER RT_kmax
 
       COMMON/radtrans_params_l/
+     &    RT_useOASIMrmud,
      &    RT_useMeanCosSolz,
      &    RT_useNoonSolz,
      &    RT_useSEAICE,
      &    RT_haveIce
+      LOGICAL RT_useOASIMrmud
       LOGICAL RT_useMeanCosSolz
       LOGICAL RT_useNoonSolz
       LOGICAL RT_useSEAICE

pkg/radtrans/radtrans_calc.F

@@ -20,6 +20,10 @@ SUBROUTINE RADTRANS_CALC(
 #include "GRID.h"
 #include "EEPARAMS.h"
 #include "PARAMS.h"
+#ifdef ALLOW_OASIM
+#include "OASIM_SIZE.h"
+#include "OASIM_FIELDS.h"
+#endif
 #include "RADTRANS_SIZE.h"
 #include "RADTRANS_PARAMS.h"
 #include "RADTRANS_FIELDS.h"
@@ -67,6 +71,9 @@ SUBROUTINE RADTRANS_CALC(
       _RL x_k(Nr), y_k(Nr)
       _RL r1_k(Nr), r2_k(Nr)
       _RL kappa1_k(Nr), kappa2_k(Nr)
+#ifdef ALLOW_OASIM
+      INTEGER loa
+#endif
 
       _RL Ed(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
       _RL Es(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
@@ -106,15 +113,18 @@ SUBROUTINE RADTRANS_CALC(
 C ======================================================================
 C--   compute solar zenith angles
 
-      IF (RT_useMeanCosSolz) THEN
-C     from day-time average of cosine of zenith angle
+C     OASIM computes its own
+      IF (.NOT.(useOASIM .AND. RT_useOASIMrmud)) THEN
+
+       IF (RT_useMeanCosSolz) THEN
+C      from day-time average of cosine of zenith angle
 
         CALL RADTRANS_DECLINATION_SPENCER(delta, myTime, myIter, myThid)
         CALL RADTRANS_SOLZ_DAYTIME(solz, delta, bi, bj,
      &                             iMin, iMax, jMin, jMax, myThid)
 
-      ELSE
-C     compute zenith angle at local noon
+       ELSE
+C      compute zenith angle at local noon
 
         IF (RT_useNoonSolz) THEN
          isec = 12*3600
@@ -128,10 +138,13 @@ SUBROUTINE RADTRANS_CALC(
      I                    myTime, myIter, myThid )
 #endif
 
-      ENDIF
+       ENDIF
 
-c     compute 1/cos(zenith) for direct light below surface
-      CALL RADTRANS_RMUD_BELOW(rmud,solz,iMin,iMax,jMin,jMax,myThid)
+c      compute 1/cos(zenith) for direct light below surface
+       CALL RADTRANS_RMUD_BELOW(rmud,solz,iMin,iMax,jMin,jMax,myThid)
+
+C     not useOASIM
+      ENDIF
 
 C ======================================================================
       DO j=jMin,jMax
@@ -145,14 +158,41 @@ SUBROUTINE RADTRANS_CALC(
 
         DO l = 1,nlam
 
+#ifdef ALLOW_OASIM
+         IF (useOASIM) THEN
+
+C use light computed in the oasim package
+          DO loa = 1,nlt
+           IF (nint(RT_wbRefWLs(l)) .EQ. oasim_lam(loa)) EXIT
+          ENDDO
+          IF (loa .gt. nlt) THEN
+            print*,'RADTRANS_CALC: waveband not found in OASIM:',
+     &             RT_wbRefWLs(l)
+            STOP 'RADTRANS_CALC: waveband not found in OASIM'
+          ENDIF
+          Edwsf = oasim_edbelow(i,j,bi,bj,loa)*RT_oasimWgt(l)
+          Eswsf = oasim_esbelow(i,j,bi,bj,loa)*RT_oasimWgt(l)
+
+          IF (RT_useOASIMrmud) THEN
+           rmud(i,j) = OASIM_rmud(i,j,bi,bj)
+          ENDIF
+
+         ELSE
+#else
+         IF (.TRUE.) THEN
+#endif /* ALLOW_OASIM */
+
 C use read-in light
           Edwsf = RT_Ed_sfc(i,j,bi,bj,l)
           Eswsf = RT_Es_sfc(i,j,bi,bj,l)
 
-          IF (myiter .ge. 0) THEN
-           Edwsf = Edwsf*(1.0 _d 0 - RT_iceFrac(i,j,bi,bj))
-           Eswsf = Eswsf*(1.0 _d 0 - RT_iceFrac(i,j,bi,bj))
-          ENDIF
+C endif useOASIM
+         ENDIF
+
+         IF (myiter .ge. 0) THEN
+          Edwsf = Edwsf*(1.0 _d 0 - RT_iceFrac(i,j,bi,bj))
+          Eswsf = Eswsf*(1.0 _d 0 - RT_iceFrac(i,j,bi,bj))
+         ENDIF
 
          DO k=1,Nr
           a_k(k) = a(i,j,k,l)

pkg/radtrans/radtrans_check.F

@@ -40,7 +40,7 @@ SUBROUTINE RADTRANS_CHECK( myThid )
 
 C !LOCAL VARIABLES: ====================================================
       LOGICAL haveAny,haveAll,addBlkLn
-      INTEGER oUnit, l, iL, errCount
+      INTEGER oUnit, l, iL, errCount, cnt
       CHARACTER*(MAX_LEN_MBUF) msgBuf
       CHARACTER*1 blkLin
 
@@ -78,9 +78,9 @@ SUBROUTINE RADTRANS_CHECK( myThid )
 C     need sun package unless using mean cos(solz)
 
 #ifndef ALLOW_SUN
-      IF (.NOT.RT_useMeanCosSolz) THEN
+      IF (.NOT.(RT_useMeanCosSolz.OR.RT_useOASIMrmud)) THEN
         WRITE(msgBuf,'(2A)') 'RADTRANS_CALC: need to compile the sun ',
-     &   'package unless RT_useMeanCosSolz is used.'
+     &   'package unless RT_useMeanCosSolz or oasim pkg is used.'
         CALL PRINT_ERROR( msgBuf, myThid )
         CALL ALL_PROC_DIE( 0 )
         STOP 'ABNORMAL END: S/R RADTRANS_CHECK'
@@ -90,9 +90,28 @@ SUBROUTINE RADTRANS_CHECK( myThid )
 C ----------------------------------------------------------------------
 C     requirements for various solar zenith angle definitions
 
-      IF (RT_useNoonSolz.AND.RT_useMeanCosSolz)THEN
+      IF (RT_useOASIMrmud .AND. .NOT.useOASIM) THEN
+        WRITE(msgBuf,'(2A)') 'RADTRANS_CHECK: ',
+     &       'RT_useOASIMrmud requires the OASIM package'
+        CALL PRINT_ERROR( msgBuf, myThid )
+        CALL ALL_PROC_DIE( 0 )
+        STOP 'ABNORMAL END: S/R RADTRANS_CHECK'
+      ENDIF
+
+      IF (useOASIM .AND. .NOT.RT_useOASIMrmud) THEN
         WRITE(msgBuf,'(2A)') '** WARNING ** RADTRANS_CHECK: ',
-     &       'Only one of RT_useMeanCosSolz, RT_useNoonSolz can be set'
+     &       'RT_useOASIMrmud=.FALSE. but is recommended with useOASIM'
+        CALL PRINT_MESSAGE(msgBuf,errorMessageUnit,SQUEEZE_RIGHT,myThid)
+      ENDIF
+
+C     make sure only on solz option is set
+      cnt = 0
+      IF (RT_useOASIMrmud) cnt = cnt + 1
+      IF (RT_useMeanCosSolz) cnt = cnt + 1
+      IF (RT_useNoonSolz) cnt = cnt + 1
+      IF (cnt .GT. 1) THEN
+        WRITE(msgBuf,'(2A)') 'RADTRANS_CHECK: only one of ',
+     &  'RT_useOASIMrmud, RT_useMeanCosSolz, RT_useNoonSolz can be set'
         CALL PRINT_ERROR( msgBuf, myThid )
         CALL ALL_PROC_DIE( 0 )
         STOP 'ABNORMAL END: S/R RADTRANS_CHECK'
@@ -297,6 +316,9 @@ SUBROUTINE RADTRANS_CHECK( myThid )
       CALL WRITE_0D_I(RT_kmax,INDEX_NONE,
      &'RT_kmax =',
      &' /* maximum depth index for radtrans computations */')
+      CALL WRITE_0D_L(RT_useOASIMrmud,INDEX_NONE,
+     &'RT_useOASIMrmud =',
+     &' /* use mean solar zenith angle from OASIM */')
       CALL WRITE_0D_L(RT_useMeanCosSolz,INDEX_NONE,
      &'RT_useMeanCosSolz =',
      &' /* flag for using mean daytime cosine of solar zenith angle */')
@@ -306,6 +328,9 @@ SUBROUTINE RADTRANS_CHECK( myThid )
       CALL WRITE_0D_RL(RT_sfcIrrThresh,INDEX_NONE,
      &'RT_sfcIrrThresh =',
      &' /* minimum irradiance for radiative transfer computations */')
+      CALL WRITE_1D_RL(RT_oasimWgt,nlam,INDEX_NONE,
+     &'RT_oasimWgt =',
+     &' /* weight multiplying each waveband from OASIM */')
       CALL WRITE_0D_RL(RT_rmus,INDEX_NONE,
      &'RT_rmus =',
      &' /* mean inv. cosine of zenith angle for dwnwrd diff. irr. */')

pkg/radtrans/radtrans_readparms.F

@@ -104,9 +104,11 @@ SUBROUTINE RADTRANS_READPARMS( myThid )
      &    RT_wbEdges,
      &    RT_wbRefWLs,
      &    RT_kmax,
+     &    RT_useOASIMrmud,
      &    RT_useMeanCosSolz,
      &    RT_useNoonSolz,
-     &    RT_sfcIrrThresh
+     &    RT_sfcIrrThresh,
+     &    RT_oasimWgt
 
       NAMELIST/RADTRANS_DEPENDENT/
      &    RT_wbWidths
@@ -136,6 +138,7 @@ SUBROUTINE RADTRANS_READPARMS( myThid )
        RT_Es_exfremo_intercept(l) = 0.0 _d 0
        RT_Ed_exfremo_slope(l)     = 0.0 _d 0
        RT_Es_exfremo_slope(l)     = 0.0 _d 0
+       RT_oasimWgt(l)             = 1.0 _d 0
       ENDDO
       RT_E_period     = 0.0 _d 0
       RT_E_RepCycle   = repeatPeriod
@@ -192,6 +195,7 @@ SUBROUTINE RADTRANS_READPARMS( myThid )
         RT_wbRefWLs(l) = -1 _d 0
       ENDDO
       RT_kmax = Nr
+      RT_useOASIMrmud = useOASIM
       RT_useMeanCosSolz = .FALSE.
       RT_useNoonSolz = .FALSE.
       RT_sfcIrrThresh = 1 _d -4