Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:37:30 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
09a6f3668a Jeff*0001 #include "ctrparam.h"
                0002 #include "ATM2D_OPTIONS.h"
                0003 
                0004 C     !INTERFACE:
                0005       SUBROUTINE INIT_ATM2D(dtatm, dtocn, dtcouple, myThid )
                0006 C     *==========================================================*
                0007 C     |  INIT_1DTO2D                                             |
                0008 C     |    This initialization routine should be run after the   |
                0009 c     |    the ocean grid/pickup have been read in.               |
                0010 c     |                                                          |
                0011 c     |  Note: grid variable indices bi,bj are hard-coded 1,1    |
                0012 c     |  This should work if coupler or atmos/coupler on one     |
                0013 c     |  machine.                                                |
                0014 c     |                                                          |
                0015 C     *==========================================================*
                0016 c
                0017         IMPLICIT NONE
                0018 
                0019 C     === Global Atmosphere Variables ===
                0020 #include "ATMSIZE.h"
                0021 #include "AGRID.h"
                0022 
                0023 C     === Global Ocean Variables ===
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
                0028 
                0029 C     === Global SeaIce Parameters ===
                0030 #include "THSICE_PARAMS.h"
                0031 
                0032 C     === Atmos/Ocean/Seaice Interface Variables ===
                0033 #include "ATM2D_VARS.h"
9274434acc Jean*0034 
09a6f3668a Jeff*0035 
                0036 C     !INPUT/OUTPUT PARAMETERS:
                0037 C     === Routine arguments ===
                0038 C     dtatm, dtocn, dtcouple - Timesteps from couple.nml (hours)
                0039 C     myThid - Thread no. that called this routine.
                0040       INTEGER dtatm, dtocn, dtcouple
                0041       INTEGER myThid
                0042 
                0043 C     LOCAL VARIABLES:
                0044       INTEGER i,j,jj
0efd285817 Jeff*0045       INTEGER ib, ibj1, ibj2  ! runoff band loop counters
09a6f3668a Jeff*0046       INTEGER j_atm, mn
                0047       INTEGER dUnit
0efd285817 Jeff*0048       _RL end1, end2, enda1, enda2, enda3 !used to compute grid conv areas
                0049       _RL totrun_b(sNy) ! total file "runoff" in runoff bands
09a6f3668a Jeff*0050       _RL a1,a2
0efd285817 Jeff*0051       _RS atm_dyG(jm0)  ! southern point/(boundary) of atmos grid
                0052       DATA atm_dyG/2.0,44*4.0,2.0/  ! grid spacing for atmosphere
09a6f3668a Jeff*0053 
                0054       dtatmo = dtatm * 3600.
                0055       dtocno = dtocn * 3600.
                0056       dtcouplo= dtcouple * 3600.
                0057 
                0058 C override data.ice seaice time step parms
                0059 C these will need to change if coupling procedure changed
                0060       thSice_deltaT = dtcouplo
                0061       thsIce_dtTemp = dtatmo
                0062       ocean_deltaT = dtcouplo
                0063 
                0064 CJRS  This next check - only kill it if not MPI?
                0065       IF (dtocno.NE.dTtracerLev(1)) THEN
                0066         PRINT *,'Ocean tracer timestep differs between coupler '
                0067         PRINT *,'and the ocean data file'
                0068         STOP
                0069       ENDIF
                0070 
                0071 c Assuming the atmospheric grid array not passed, do this:
                0072       atm_yG(1)=-90.0
                0073       DO j_atm=2,jm0
                0074         atm_yG(j_atm)=atm_yG(j_atm-1)+atm_dyG(j_atm-1)
                0075         atm_yC(j_atm-1)=(atm_yG(j_atm-1)+atm_yG(j_atm))/2.0
                0076       ENDDO
                0077       atm_yC(jm0)=atm_yG(jm0)+atm_dyG(jm0)/2.0
                0078 
                0079 c end atmos grid initialization
                0080 
                0081       atm_oc_ind(1)=2
                0082       atm_oc_wgt(1)=1. _d 0
9274434acc Jean*0083       atm_oc_frac1(1)= (sin(yG(1,2,1,1)*deg2rad) -
0efd285817 Jeff*0084      &        sin(yG(1,1,1,1)*deg2rad))/
                0085      &        (sin(atm_yG(3)*deg2rad)-sin(atm_yG(1)*deg2rad))
                0086       atm_oc_frac2(1)= 0. _d 0   ! assumes ocean(1) fits in atm(1)
09a6f3668a Jeff*0087       atm_oc_ind(sNy)=jm0-1
                0088       atm_oc_wgt(sNy)=1. _d 0
9274434acc Jean*0089       atm_oc_frac1(sNy)= (sin((yG(1,sNy,1,1) +
0efd285817 Jeff*0090      &      dyG(1,sNy,1,1)/6.37D6/deg2rad)*deg2rad)-
                0091      &      sin(yG(1,sNy,1,1)*deg2rad))/
9274434acc Jean*0092      &      (sin((atm_yG(jm0)+atm_dyG(jm0))*deg2rad)-
0efd285817 Jeff*0093      &      sin(atm_yG(jm0-1)*deg2rad))
                0094       atm_oc_frac2(sNy)= 0. _d 0   ! assumes ocean(1) fits in atm(1)
                0095 
                0096       endwgt1 = sin(atm_yG(2)*deg2rad)            !hard-coded that the atmos
                0097       endwgt2 = sin(atm_yG(3)*deg2rad) - endwgt1  !grid is same in NH and SH
                0098       endwgt1 = endwgt1 + 1. _d 0                 !and goes 90S to 90N
                0099       rsumwgt = 1. _d 0/(endwgt1 + endwgt2)
                0100 
                0101       atm_yG(2)=atm_yG(1) ! grid now combined atm end points
                0102       atm_yG(jm0)=90. _d 0
09a6f3668a Jeff*0103 
                0104       DO j=2, sNy-1
                0105 
                0106         DO jj=2,jm0-1
0efd285817 Jeff*0107           IF  ((yG(1,j,1,1).GE.atm_yG(jj)).AND.
                0108      &         (yG(1,j,1,1).LT.atm_yG(jj+1))) j_atm=jj
9274434acc Jean*0109         ENDDO
                0110 
09a6f3668a Jeff*0111         atm_oc_ind(j)=j_atm
0efd285817 Jeff*0112         end1= sin(yG(1,j,1,1) *deg2rad)
                0113         end2= sin(yG(1,j+1,1,1) *deg2rad)
                0114         enda1 = sin(atm_yG(j_atm) *deg2rad)
                0115         enda2 = sin(atm_yG(j_atm+1) *deg2rad)
                0116         IF ( yG(1,j+1,1,1) .GT. atm_yG(j_atm+1) ) THEN
                0117            enda3 = sin(atm_yG(j_atm+2) *deg2rad)
                0118            atm_oc_wgt(j)=(enda2-end1)/ (end2-end1)
                0119            atm_oc_frac1(j)= (enda2-end1) / (enda2 - enda1)
                0120            atm_oc_frac2(j)= (end2 - enda2) / (enda3 - enda2)
09a6f3668a Jeff*0121         ELSE
                0122           atm_oc_wgt(j)=1. _d 0
0efd285817 Jeff*0123           atm_oc_frac1(j)= (end2-end1)/ (enda2-enda1)
                0124           atm_oc_frac2(j)=0. _d 0
09a6f3668a Jeff*0125         ENDIF
                0126       ENDDO
f67c0bf3c1 Jeff*0127 
                0128 C     compute tauv interpolation points
                0129       tauv_jpt(1) = 2         ! south pole point; s/b land
                0130       tauv_jwght(1) = 1. _d 0
                0131       DO j=2, sNy
                0132         DO jj=1,jm0-1
                0133           IF (( yG(1,j,1,1) .GE. atm_yC(jj)).AND.
                0134      &        ( yG(1,j,1,1) .LT. atm_yC(jj+1))) j_atm=jj
                0135         ENDDO
                0136         tauv_jpt(j)= j_atm
7b7831d041 Jean*0137         tauv_jwght(j)= 1. _d 0 - (yG(1,j,1,1) - atm_yC(j_atm)) /
f67c0bf3c1 Jeff*0138      &                 (atm_yC(j_atm+1) - atm_yC(j_atm))
                0139       ENDDO
                0140 
                0141 C      DO j=1,sNy
                0142 C      print *, 'j, tauv_jpt:', j,tauv_jpt(j),tauv_jwght(j)
                0143 C      ENDDO
                0144 
09a6f3668a Jeff*0145 c
                0146 c find land fraction
                0147 c
                0148       DO j_atm=1,jm0
                0149         cflan(j_atm)=0. _d 0
                0150         ocnArea(j_atm)=0. _d 0
                0151       ENDDO
                0152 
                0153       DO j=1,sNy
                0154         DO i=1,sNx
                0155           IF (maskC(i,j,1,1,1).EQ.1.) THEN
                0156             ocnArea(atm_oc_ind(j))=ocnArea(atm_oc_ind(j)) +
                0157      &                           rA(i,j,1,1)*atm_oc_wgt(j)
0efd285817 Jeff*0158             IF (atm_oc_wgt(j).LT.1.d0) THEN
09a6f3668a Jeff*0159               ocnArea(atm_oc_ind(j)+1)=ocnArea(atm_oc_ind(j)+1) +
                0160      &                           rA(i,j,1,1)*(1.d0-atm_oc_wgt(j))
                0161             ENDIF
                0162           ENDIF
                0163         ENDDO
                0164       ENDDO
                0165 
                0166       DO j_atm=3,jm0-2
9274434acc Jean*0167         cflan(j_atm)=1. _d 0 - ocnArea(j_atm) /
                0168      &              (2. _d 0 * PI * 6.37 _d 6 * 6.37 _d 6 *
09a6f3668a Jeff*0169      &    (sin(atm_yG(j_atm+1)*deg2rad) - sin(atm_yG(j_atm)*deg2rad)))
                0170         if (cflan(j_atm).LT.1. _d -14) cflan(j_atm)=0. _d 0
                0171       ENDDO
                0172 
                0173 C     deal with the combined atmos grid end cells...
9274434acc Jean*0174       cflan(2)= 1. _d 0 - ocnArea(2) /
09a6f3668a Jeff*0175      &         (2. _d 0*PI*6.37 _d 6*6.37 _d 6*
                0176      &         (sin(atm_yG(3)*deg2rad)+1. _d 0))
                0177       IF (cflan(2).LT.1. _d -14) cflan(2)=0. _d 0
                0178       cflan(1)=cflan(2)
9274434acc Jean*0179       cflan(jm0-1)= 1.d0 - ocnArea(jm0-1) /
09a6f3668a Jeff*0180      &             (2. _d 0*PI*6.37 _d 6*6.37 _d 6*
                0181      &             (1. _d 0-sin(atm_yG(jm0-1)*deg2rad)))
                0182       IF (cflan(jm0-1).LT.1. _d -14) cflan(jm0-1)=0. _d 0
                0183       cflan(jm0)=cflan(jm0-1)
                0184 
9274434acc Jean*0185       PRINT *,'Land fractions on atmospheric grid: '
09a6f3668a Jeff*0186       PRINT *, cflan
                0187       PRINT *,'Lookup grid index, weights:'
                0188       PRINT *, atm_oc_ind,atm_oc_wgt
0efd285817 Jeff*0189 C      PRINT *,'Lookup fraction 1 of atmos grid:'
                0190 C      PRINT *, atm_oc_frac1
                0191 C      PRINT *,'Lookup fraction 2 of atmos grid:'
                0192 C      PRINT *, atm_oc_frac2
09a6f3668a Jeff*0193 
                0194 c
                0195 c read in mean 1D atmos wind files -- store in memory
                0196 c
                0197       DO j_atm=1,jm0
                0198         DO mn=1,nForcingPer
                0199           atau(j_atm,mn)=0. _d 0
                0200           atav(j_atm,mn)=0. _d 0
                0201           awind(j_atm,mn)=0. _d 0
                0202         ENDDO
                0203       ENDDO
                0204 
                0205       CALL MDSFINDUNIT( dUnit, myThid )
                0206 
                0207       IF ( atmosTauuFile .NE. ' '  ) THEN
                0208          OPEN(dUnit, FILE=atmosTauuFile,STATUS='old',
9274434acc Jean*0209      &        ACCESS='direct', RECL=8*jm0*nForcingPer,
09a6f3668a Jeff*0210      &        FORM='unformatted')
                0211          READ(dUnit,REC=1), atau
                0212          CLOSE(dUnit)
                0213       ENDIF
                0214 
                0215       IF ( atmosTauvFile .NE. ' '  ) THEN
                0216          OPEN(dUnit, FILE=atmosTauvFile, STATUS='old',
9274434acc Jean*0217      &        ACCESS='direct', RECL=8*jm0*nForcingPer,
09a6f3668a Jeff*0218      &        FORM='unformatted')
                0219          READ(dUnit, REC=1), atav
                0220          CLOSE(dUnit)
                0221       ENDIF
                0222 
                0223       IF ( atmosWindFile .NE. ' '  ) THEN
                0224          OPEN(dUnit, FILE=atmosWindFile, STATUS='old',
9274434acc Jean*0225      &        ACCESS='direct', RECL=8*jm0*nForcingPer,
09a6f3668a Jeff*0226      &        FORM='unformatted')
                0227          READ(dUnit, REC=1), awind
                0228          CLOSE(dUnit)
                0229       ENDIF
                0230 
                0231 C The polar data point values for winds are effectively N/A given the
                0232 C pole issue... although they are read in here, they are never used in
                0233 C any calculations, as the polar ocean points access the data at atmos
                0234 C 2 and jm0-1 points.
                0235 
                0236 
                0237 c read in runoff data
                0238 c to put runoff into specific grid cells
                0239 c
                0240       IF ( runoffMapFile .NE. ' ' ) THEN
                0241         CALL READ_FLD_XY_RL( runoffMapFile, ' ',
                0242      &                      runoffVal, 0, myThid )
                0243       ELSE
                0244         DO j=1,sNy
                0245           DO i=1,sNx
                0246             if ( (maskC(i,j,1,1,1).EQ.1.) .AND.
                0247      &          ( (maskC(i-1,j,1,1,1).EQ.0.).OR.
                0248      &            (maskC(i+1,j,1,1,1).EQ.0.).OR.
                0249      &            (maskC(i,j-1,1,1,1).EQ.0.).OR.
                0250      &            (maskC(i,j+1,1,1,1).EQ.0.).OR.
                0251      &            (maskC(i+1,j+1,1,1,1).EQ.0.).OR.
                0252      &            (maskC(i-1,j-1,1,1,1).EQ.0.).OR.
                0253      &            (maskC(i+1,j-1,1,1,1).EQ.0.).OR.
                0254      &            (maskC(i-1,j+1,1,1,1).EQ.0.) ) ) THEN
                0255               runoffVal(i,j)=1. _d 0
                0256             ENDIF
                0257           ENDDO
                0258         ENDDO
                0259       ENDIF
                0260 
                0261       DO ib=1,numBands
                0262         ibj1=1
                0263         if (ib.GT.1) ibj1=  rband(ib-1)+1
                0264         ibj2=sNy
                0265         if (ib.LT.numBands) ibj2= rband(ib)
                0266         totrun_b(ib)=0.d0
                0267         DO j=ibj1,ibj2
                0268           DO i=1,sNx
                0269             totrun_b(ib)=totrun_b(ib)+runoffVal(i,j)*maskC(i,j,1,1,1)
                0270           ENDDO
                0271         ENDDO
                0272         DO j=ibj1,ibj2
                0273           runIndex(j)= ib     ! for lookup of rband as fn. of latitude
                0274           DO i=1,sNx
                0275             runoffVal(i,j)=runoffVal(i,j)*maskC(i,j,1,1,1)/totrun_b(ib)
                0276           ENDDO
                0277         ENDDO
                0278       ENDDO
                0279 
                0280       CALL INIT_SUMVARS(myThid)
9274434acc Jean*0281 
                0282 C     Initialize 1D diagnostic variables
09a6f3668a Jeff*0283       DO j_atm=1,jm0
                0284         DO mn=1,nForcingPer
                0285           sum_tauu_ta(j_atm,mn)= 0. _d 0
                0286           sum_tauv_ta(j_atm,mn)= 0. _d 0
                0287           sum_wsocean_ta(j_atm,mn)= 0. _d 0
                0288           sum_ps4ocean_ta(j_atm,mn)= 0. _d 0
                0289         ENDDO
                0290       ENDDO
                0291 
9274434acc Jean*0292 C     Initialize 2D diagnostic variables
09a6f3668a Jeff*0293       DO i=1-OLx,sNx+OLx
                0294         DO j=1-OLy,sNy+OLy
                0295           DO mn=1,nForcingPer
                0296             qnet_atm_ta(i,j,mn)= 0. _d 0
                0297             evap_atm_ta(i,j,mn)= 0. _d 0
                0298             precip_atm_ta(i,j,mn)= 0. _d 0
                0299             runoff_atm_ta(i,j,mn)= 0. _d 0
                0300             sum_qrel_ta(i,j,mn)= 0. _d 0
                0301             sum_frel_ta(i,j,mn)= 0. _d 0
                0302             sum_iceMask_ta(i,j,mn)= 0. _d 0
                0303             sum_iceHeight_ta(i,j,mn)= 0. _d 0
                0304             sum_iceTime_ta(i,j,mn)= 0. _d 0
                0305             sum_oceMxLT_ta(i,j,mn)= 0. _d 0
                0306             sum_oceMxLS_ta(i,j,mn)= 0. _d 0
                0307           ENDDO
                0308           qnet_atm(i,j)= 0. _d 0
                0309           evap_atm(i,j)= 0. _d 0
                0310           precip_atm(i,j)= 0. _d 0
                0311           runoff_atm(i,j)= 0. _d 0
                0312           sum_qrel(i,j)= 0. _d 0
                0313           sum_frel(i,j)= 0. _d 0
                0314           sum_iceMask(i,j)= 0. _d 0
                0315           sum_iceHeight(i,j)= 0. _d 0
                0316           sum_iceTime(i,j)= 0. _d 0
                0317           sum_oceMxLT(i,j)= 0. _d 0
                0318           sum_oceMxLS(i,j)= 0. _d 0
                0319         ENDDO
                0320       ENDDO
                0321 
8e101fde6e Jeff*0322 C     Initialize year-end diags and max/min seaice variables
                0323       SHice_min = 1. _d 18
                0324       NHice_min = 1. _d 18
                0325       SHice_max = 0. _d 0
                0326       NHice_max = 0. _d 0
                0327       sst_tave=  0. _d 0
                0328       sss_tave=  0. _d 0
                0329       HF2ocn_tave=  0. _d 0
                0330       FW2ocn_tave=  0. _d 0
678e639be8 Jeff*0331       CO2flx_tave=  0. _d 0
8e101fde6e Jeff*0332       OPEN(25,FILE='resocean.dat',STATUS='replace')
                0333       CLOSE(25)
                0334 
09a6f3668a Jeff*0335 C     Initialize following for safety and/or cold start
                0336       DO i=1-OLx,sNx+OLx
                0337         DO j=1-OLy,sNy+OLy
                0338           pass_runoff(i,j)= 0. _d 0
                0339           pass_qnet(i,j)= 0. _d 0
                0340           pass_evap(i,j)= 0. _d 0
                0341           pass_precip(i,j)= 0. _d 0
                0342           pass_fu(i,j)= 0. _d 0
                0343           pass_fv(i,j)= 0. _d 0
                0344           pass_wspeed(i,j)= 0. _d 0
                0345           pass_solarnet(i,j)= 0. _d 0
                0346           pass_slp(i,j)= 0. _d 0
                0347           pass_siceLoad(i,j)= 0. _d 0
                0348           pass_pCO2(i,j)= 0. _d 0
7b7831d041 Jean*0349           pass_prcAtm(i,j) = 0. _d 0
                0350           snowPrc    (i,j) = 0. _d 0
09a6f3668a Jeff*0351           sFluxFromIce(i,j)= 0. _d 0
                0352         ENDDO
                0353       ENDDO
                0354 
678e639be8 Jeff*0355 C     Initialize following (if ocn carbon not passed)
                0356       DO i=1,sNx
                0357         DO j=1,sNy
                0358           oFluxCO2(i,j) = 0. _d 0
                0359         ENDDO
                0360       ENDDO
                0361 
09a6f3668a Jeff*0362       RETURN
                0363       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!