C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/AITCZ/code/mitphys_do_atmos_physics.F,v 1.1 2003/08/20 15:24:59 czaja Exp $
C     $Name:  $

#include "CPP_OPTIONS.h"
#undef OLD_AIM_GRIG_MAPPING

      SUBROUTINE MITPHYS_DO_ATMOS_PHYSICS( phi_hyd, currentTime, myThid)
C     /==================================================================\
C     | S/R MITPHYS_DO_ATMOS_PHYSICS                                     |
C     |==================================================================|
C     | Interface interface between atmospheric physics package and the  |
C     | dynamical model.                                                 |
C     | Routine calls physics pacakge after mapping model variables to   |
C     | the package grid. Package should derive and set tendency terms   |
C     | which can be included as external forcing terms in the dynamical |
C     | tendency routines. Packages should communicate this information  |
C     | through common blocks.                                           |
C     \==================================================================/
      IMPLICIT NONE

C     -------------- Global variables ------------------------------------
C     Physics package
#include "atparam.h"
#include "atparam1.h"
#include "MITPHYS_PARAMS.h"
      INTEGER NGP
      INTEGER NLON
      INTEGER NLAT
      INTEGER NLEV
      PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
#include "thermo.h"
#include "com_mitphysvar.h"
#include "com_forcing1.h"
#include "Lev_def.h"
C     MITgcm
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"

C     -------------- Routine arguments -----------------------------------
      _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL currentTime
      INTEGER myThid

C omp modif: shapiro filter on the convective mass flux:
C      _RL cbmf_tmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL dist


#ifdef ALLOW_MITPHYS
C     -------------- Local variables -------------------------------------
C     I,J,K,I2,J2   - Loop counters
C     tYear         - Fraction into year
C     mnthIndex     - Current month
C     prevMnthIndex - Month last time this routine was called.
C     tmp4          - I/O buffer ( 32-bit precision )
C     fNam          - Work space for file names
C     mnthNam       - Month strings
C     hInital       - Initial height of pressure surfaces (m)
C     pSurfs        - Pressure surfaces (Pa)
C     Katm          - Atmospheric K index
      INTEGER I
      INTEGER I2
      INTEGER J
      INTEGER J2
      INTEGER K
      INTEGER IG0
      INTEGER JG0
      _RL    tYear
      INTEGER mnthIndex
      INTEGER prevMnthIndex
      DATA    prevMnthIndex / 0 /
      SAVE    prevMnthIndex
      LOGICAL FirstCall
      DATA    FirstCall /.TRUE./
      SAVE    FirstCall
      LOGICAL CALLFirst
      DATA    CALLFirst /.TRUE./
      SAVE    CALLFirst
      INTEGER nxIo 
      INTEGER nyIo
      PARAMETER ( nxIo = 128, nyIo = 64 )
      _RL  tmp4(nxIo,nyIo)
      CHARACTER*16 fNam
      CHARACTER*3 mnthNam(12)
      DATA mnthNam /
     & 'jan', 'feb', 'mar', 'apr', 'may', 'jun', 
     & 'jul', 'aug', 'sep', 'oct', 'nov', 'dec' /
      SAVE mnthNam
      _RL    hInitial(Nr)
      _RL    hInitialW(Nr)
      _RL    pSurfs(Nr)
c 5 levels      DATA    pSurfs   / 950.D2,775.D2, 500.D2,  250.D2, 75.D2 /
      DATA    pSurfs  / 1000.0D2, 975.0D2, 950.0D2, 925.0D2, 900.0D2,
     :                   875.0D2, 850.0D2, 825.0D2, 800.0D2, 775.0D2,
     :                   750.0D2, 725.0D2, 700.0D2, 675.0D2, 650.0D2,
     :                   625.0D2, 600.0D2, 575.0D2, 550.0D2, 525.0D2,
     :                   500.0D2, 475.0D2, 450.0D2, 425.0D2, 400.0D2,
     :                   375.0D2, 350.0D2, 325.0D2, 300.0D2, 275.0D2,
     :                   250.0D2, 225.0D2, 200.0D2, 175.0D2, 150.0D2,
     :                   125.0D2, 100.0D2,  75.0D2,  50.0D2,  25.0D2 /
      SAVE    pSurfs
      _RL Soilqmax
      _RL pGround
      INTEGER bi, bj
      INTEGER Katm

      INTEGER NEXT_CALL
      DATA NEXT_CALL /0/
      SAVE NEXT_CALL


      _RL LAT_CYCLE

C      LOGICAL DO_OMP_SHAP
C

CC (acz) check for SST
CC         write(6,*) 'beginning1 of do_atmos_phys: SST_REF=',SST_REF(1)

      bi  = 1
      bj  = 1

      IF (CurrentTime .LT. 10.d0) THEN
         DO J=1,sNy 
            DO I=1,sNx 
               I2 = (sNx)*(J-1)+I
               SST1 (I2)= SST_REF(I2)
               STL1 (I2)= STL_REF(I2)
               CBMF_DYN(I,J,bi,bj) = CBMFG1(I2)
               SST_DYN(I,J,bi,bj) = SST_REF(I2)
               STL_DYN(I,J,bi,bj) = STL_REF(I2)
            END DO
         END DO
      END IF

CC (acz) check for SST
CC         write(6,*) 'beginning2 of do_atmos_phys: SST1=',SST1(1)
CC         write(6,*) 'beginning2 of do_atmos_phys: SST_DYN=',SST_DYN(1,1,bi,bj)

      IF (NEXT_CALL .EQ. 0 ) THEN

      pGround = 1012.5D2

C     Assume only one tile per proc. for now
      bi  = 1
      bj  = 1
      IG0 = myXGlobalLo
      JG0 = myYGlobalLo

C      
C     Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.
C     Internal index mapping is linear in X and Y with a second
C     dimension for the vertical.

C     Adjustment for heave due to mean heating/cooling
C     ( I don't think the old formula was strictly "correct" for orography
C       but I have implemented it as was for now. As implemented
C       the mean heave of the bottom (K=Nr) level is calculated rather than
C       the mean heave of the base of the atmosphere. )

c -- sb: remove this with the MIT physics:
c sb c_jmc: Because AIM physics LSC is not applied in the stratosphere (top level),
c sb c      ==> move water wapor from the stratos to the surface level.
c sb       DO J = 1-Oly, sNy+Oly
c sb        DO I = 1-Olx, sNx+Olx
c sb c       k = k_surf(i,j,bi,bj)
c sb c       salt(I,J,k,bi,bj) = salt(I,J,k,bi,bj)
c sb c    &   + maskC(i,j,Nr,bi,bj)*salt(I,J,Nr,bi,bj)*drF(Nr)*recip_drF(k)
c sb        salt(I,J,Nr,bi,bj) = 0.
c sb       ENDDO
c sb      ENDDO
c sb --

C     Note the mapping here is only valid for one tile per proc.
      DO J = 1, sNy
        DO I = 1, sNx
          I2 = (sNx)*(J-1)+I
C omp: done in mitphys_do_inphys now
C latitude and logitude, assuming a spherical coordinate system
C Must be modified for cartesian coordinates
          LAT_G(I2) = yC(I,J,bi,bj)
          LON_G(I2) = xC(I,J,bi,bj)
C omp: new addition for restart files.
        CBMFG1(I2) = CBMF_DYN (I,J,bi, bj)
        SST1(I2) = SST_DYN (I,J,bi, bj)
        STL1(I2) = STL_DYN (I,J,bi, bj)

CC (acz) check for SST
CC       IF (I2 .EQ. 1) THEN
CC        write(6,*) 'SST1 set to SST_DYN at the begin. of do_atmos_phys: SST1,SSTDYN', SST1(I2),SST_DYN (1,1,bi, bj)
CC       END IF


          DO K = 1, Nr
             Katm = K 

C         UG1(I2,Katm)   = 0.5*(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))
C         VG1(I2,Katm)   = 0.5*(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj))

C         UG1(I2,Katm)   = uVel(I,J,K,bi,bj)
C         VG1(I2,Katm)   = vVel(I,J,K,bi,bj)

             UGW(I2,Katm)   = uVel(I,J,K,bi,bj)
             UGE(I2,Katm)   = uVel(I+1,J,K,bi,bj)
             VGS(I2,Katm)   = vVel(I,J,K,bi,bj)
             VGN(I2,Katm)   = vVel(I,J+1,K,bi,bj)

             WG1(I2, Katm) = wVel (I,J,K,bi, bj)

C        Physics works with temperature - not potential temp.
             TG1(I2,Katm)   = theta(I,J,K,bi,bj)/
     &            ((pGround/pSurfs(K))**(RD/CPD))
c_jmc    QG1(I2,Katm)   = salt(I,J,K,bi,bj)
             QG1(I2,Katm)   = MAX(salt(I,J,K,bi,bj), 0. _d 0)
             PHIG1(I2,Katm) = phi_hyd(I,J,K,bi,bj) + etaN(I,J,bi,bj) ! not used in MIT physics

        ENDDO

       ENDDO
      ENDDO

C      
      DO J=1,sNy
       DO I=1,sNx
        I2 = (sNx)*(J-1)+I
         PSG1(I2) = pGround
       ENDDO
      ENDDO
C      write(*,*) 'CBMF_DYN before:', CBMF_DYN 
C      write(*,*) 'CBMFGE before:', CBMFG1 

C
C
C     Physics package needs to know time of year as a fraction

C not used anymore -omp

      tYear = currentTime/(86400.*360.) -
     &        FLOAT(INT(currentTime/(86400.*360.)))
C_EqCh: Fix solar insolation with Sun directly overhead on Equator
      tYear = 0.25
C
C     Load external data needed by physics package
C     1. Albedo
C     2. Surface temperatures
C     3. Soil moisture
C     4. Snow depth            - assume no snow for now
C     5. Sea ice               - assume no sea ice for now
C     6. Land sea mask         - infer from exact zeros in soil moisture dataset
C     7. Surface geopotential  - to be done when orography is in
C                                dynamical kernel. Assume 0. for now.
      mnthIndex = INT(tYear*12.)+1
C      IF ( mnthIndex .NE. prevMnthIndex .OR.
C     &     FirstCall ) THEN
C       prevMnthIndex = mnthIndex

C 1. albedo (not used)

       DO J=1,sNy
        DO I=1,sNx
         I2 = (sNx)*(J-1)+I
         alb0(I2) = 0.

        ENDDO
       ENDDO

C modif omp: introduce an annual cycle in the prescribed SST.
C modif acz: introduce an annual cycle in the prescribed STL.

       IF (ANNUAL_CYCLE) THEN

C (acz) LAT_CYCLE = 1 at t=0 to be consistent with do_inphys.F
C This means that the integration starts in late summer

         LAT_CYCLE = 0.5 + 0.5 * cos(2 * pi *  
     &          (currentTime / 24./ 3600.) / 360.)


C 2. SST / STL
C modif omp: the background SST is defined thourgh the MITPHYS namelist
C modif acz: introduce analytical expression for STL

C omp: for spherical geometry
       IF (UsingSphericalPolarGrid) THEN

        DO J=1,sNy 
         DO I=1,sNx 
          I2 = (sNx)*(J-1)+I

          SST_REF(I2) = SST_BACK 
          STL_REF(I2) = SST_BACK 

C Equator to pole SST gradient

          SST_REF(I2) = SST_REF(I2) - DELT_EQ_PO  * 
     &         sin( LAT_G(I2) * pi  / 180. ) **2 
   

C Equator to northern boundary SST gradient - work only for a single tile

          SST_REF(I2) = SST_REF(I2) - DELT_EQ_BND  * 
     &         sin( LAT_G(I2) / abs(Phimin)) **2 

C Lindzen and Hou type gradient:

          SST_REF(I2) = SST_REF(I2) -   DELT_LH  * ( 
     &              3  *  ( sin ( PI * LAT_G(I2) / 180.d0 ) 
     &              - sin(PI * LAT_CYCLE * LAT_LH / 180.d0) 
     &                ) ** 2) 
   

CC (acz) SST used for TOSA1:
CC Zonally symmetric SST perturbation


         dist = SQRT ( (LAT_G(I2) - LAT_PERT*LAT_CYCLE) ** 2)

          if (dist .LE. DEL_LAT) then      
              SST_REF(I2) = SST_REF(I2) + DELT_PERT 
     &            * cos ( pi * dist / 2./DEL_LAT ) **2
          end if  

CC (acz) STL used for TOSA1:
CC
CC 1 - Warm ring centered at the equator with same
CC meridional extent (DEL_LAT) than the above (TOSA1) SST
CC The profile is further time-dependent: 
C  strongest at equinox, zero at solstice

         dist = ABS ( LAT_G(I2) )

          if (dist .LE. DEL_LAT) then
              STL_REF(I2) = STL_REF(I2) + DELT_STL_EQU
     &            * cos ( pi * dist / 2./DEL_LAT ) **2
     &            * cos ( 2 * pi * (currentTime / 24./ 3600.) 
     &                    / 360. + 0.5 * pi ) **2
          end if

C 2- Opposite sign variations north and south of equator
C North is warm at t=0

         dist = ABS ( LAT_G(I2) )

          if (dist .LE. DEL_LAT_STL) then
              STL_REF(I2) = STL_REF(I2) + DELT_STL_HEM
     &            * sin ( pi * LAT_G(I2) / DEL_LAT_STL ) 
     &            * cos ( 2 * pi * (currentTime / 24./ 3600.) / 360.)
          end if

         ENDDO
        ENDDO

      ELSE

C Cartesian Grid
        DO J=1,sNy 
         DO I=1,sNx 
          I2 = (sNx)*(J-1)+I
          SST_REF(I2) = SST_BACK 
          SST_REF(I2) = SST_REF(I2) - DELT_EQ_BND  * 
     &         sin( 0.5 * PI *  yC(I,J,bi,bj) / abs(phiMin) ) **2 

         END DO
        END DO

      END IF

      END IF
CC end for the ANNUAL_CYCLE end if loop

      IF (PRESC_SST) THEN
         DO I2=1,NGP 
            SST1(I2) = SST_REF(I2) 
         END DO
      END IF

      IF (PRESC_STL) THEN
        DO I2=1,NGP 
           STL1(I2) = STL_REF(I2)
        END DO
      END IF
       
CC (acz) check for SST
CC         write(6,*) 'end of do_atmos_physic: SST1=',SST1(1)


C 3. Soil moisture (not used)

       DO J=1,sNy
        DO I=1,sNx
         I2 = (sNx)*(J-1)+I
         soilq1(I2) = 0.
        ENDDO
       ENDDO

       Soilqmax=20.

C      ENDIF
C

C 4. Snow depth (not used)
      IF ( FirstCall ) THEN
C      Set snow depth, sea ice to zero for now
C      Land-sea mask ( figure this out from where 
C                      soil moisture is exactly zero ).
CC
CC(acz) VERY RISKY lines !!!!!!!!!
CC       DO J=1,sNy
CC        DO I=1,sNx
CC         I2 = (sNx)*(J-1)+I
CC         fMask1(I2) = 1.
CC         IF ( soilq1(I2) .EQ. 0. ) fMask1(I2) = 0.
CC         oice1(I2) = 0.
CC         snow1(I2) = 0.
CC        ENDDO
CC       ENDDO

      ENDIF
C
C Addition may 15 . Reset humidity to 0. if negative
C ---------------------------------------------------
Caja  DO K=1,Nr
Caja   DO J=1-OLy,sNy+OLy
Caja    DO I=1-Olx,sNx+OLx
Caja     IF ( salt(i,j,k,bi,bj) .LT. 0. .OR. K .EQ. Nr ) THEN
Caja      salt(i,j,k,bi,bj) = 0.
Caja     ENDIF
Caja    ENDDO
Caja   ENDDO
Caja  ENDDO
 
cccc      CALL PDRIVER( tYear )
CCC (acz) also known as mitphy_driver.F !!!
CCC 
      write(*,*) 'currentTime: ',currentTime
      CALL MPDRIVER( currentTime, float(PHYS_CALL) * deltaT ) ! Call the MIT physics
C
#ifdef ALLOW_TIMEAVE
C     Calculate diagnostics for MITPHYS
      CALL MITPHYS_CALC_DIAGS( bi, bj, currentTime, myThid )
#endif /* ALLOW_TIMEAVE */
C
      FirstCall = .FALSE.

C
C The dry adiabatic adjustment done in the convection scheme may
C have affected the temperature and humidity profiles:
C
      DO K = 1, Nr
       DO J = 1, sNy
c!       DO J = 1, sNy-1 ! because last row of latitude not meaningful here
        DO I = 1, sNx
         I2 = (sNx)*(J-1)+I
         Katm = K
         theta(I,J,K,bi,bj) = TG1(I2,Katm)
     :       * ( (pGround/pSurfs(K))**(RD/CPD) )
         salt(I,J,K,bi,bj) = QG1(I2,Katm)
C omp: dry adjustment on velocity is now done as a time-tendency term in external_forcing.
C         uVel(I,J,K,bi,bj) = UGW(I2,Katm)   
C         vVel(I,J,K,bi,bj) = VGS(I2,Katm) 
        ENDDO
       ENDDO
      END DO


       DO J = 1, sNy
        DO I = 1, sNx
         I2 = (sNx)*(J-1)+I
         CBMF_DYN(I,J,bi,bj) = CBMFG1(I2)
         SST_DYN(I,J,bi,bj) = SST1(I2)
         STL_DYN(I,J,bi,bj) = STL1(I2)
       END DO
      END DO

C -- Following section is neccesary since we are adjusting theta/salt
C    directly and the overlaps need to be updated accordingly.
      _EXCH_XYZ_R8(theta,myThid)
      _EXCH_XYZ_R8(salt,myThid)
      _EXCH_XY_R8(CBMF_DYN,myThid)
      _EXCH_XY_R8(SST_DYN,myThid)
      _EXCH_XY_R8(STL_DYN,myThid)


C omp: shapiro-filtering the convective mass flux:
C -still assuming single processor:
C omp: dropped... to be implemented again.

c$$$
c$$$      DO_OMP_SHAP = .FALSE.
c$$$      IF (DO_OMP_SHAP)THEN
c$$$         DO J = 1, sNy
c$$$            DO I = 1, sNx
c$$$               I2 = I+(J-1)*sNx
c$$$               cbmf_tmp(i,j,bi,bj) = CBMFG1(I2)
c$$$            end do
c$$$         end do
c$$$         
c$$$         CALL OMP_SHAP_FILTER(cbmf_tmp,currentTime, myThid)
c$$$         
c$$$         DO J = 1, sNy
c$$$            DO I = 1, sNx
c$$$               I2 = I+(J-1)*sNx
c$$$               CBMFG1(I2) = CBMF_tmp(i,j,bi,bj)
c$$$            end do
c$$$         end do
c$$$      END IF

C
C
C     For velocity tendencies on C grid need to interpolate
C     to do this in multi-processor mode we copy U and V tendencies
C     into dynamics conforming array and then exchange.
C     ** NOTE - Exchange at this point is not compatible with
C               multiple-tiles per compute thread/process.
C               However, AIM code itself is not compatible with
C               this as its global data assumes only one "tile".
c sb      CALL AIM_UTVT2DYN( bi, bj, myThid )
c sb      _EXCH_XYZ_R8( AIM_UT, myThid )
c sb      _EXCH_XYZ_R8( AIM_VT, myThid )

c sb: remove this for the moment (no drag computed by MITPHYS):
c sb      DO J = 1, sNy
c sb       DO I = 1, sNx
c sb         I2 = I+(J-1)*sNx
c sb         aim_drag(I,J,bi,bj) = DRAG(I2)
c sb       ENDDO
c sb      ENDDO
c sb      _EXCH_XY_R8( aim_drag, myThid )
C
      NEXT_CALL = PHYS_CALL
C      write(*,*) 'CBMF_DYN after:', CBMF_DYN
c      write(*,*) 'CBMFG1 after:', CBMFG1 
C       write(*,*) 'LAT :', LAT_G
C       write(*,*) 'LON :', LON_G
      END IF
C(acz) end if for the NEXT_CALL if
      NEXT_CALL = NEXT_CALL -1

#endif /* ALLOW_MITPHYS */

      RETURN
      END