C $Header: /u/gcmpack/MITgcm/pkg/atm2d/init_atm2d.F,v 1.8 2013/05/02 20:37:52 jmc Exp $ C $Name: $ #include "ctrparam.h" #include "ATM2D_OPTIONS.h" C !INTERFACE: SUBROUTINE INIT_ATM2D(dtatm, dtocn, dtcouple, myThid ) C *==========================================================* C | INIT_1DTO2D | C | This initialization routine should be run after the | c | the ocean grid/pickup have been read in. | c | | c | Note: grid variable indices bi,bj are hard-coded 1,1 | c | This should work if coupler or atmos/coupler on one | c | machine. | c | | C *==========================================================* c IMPLICIT NONE C === Global Atmosphere Variables === #include "ATMSIZE.h" #include "AGRID.h" C === Global Ocean Variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" C === Global SeaIce Parameters === #include "THSICE_PARAMS.h" C === Atmos/Ocean/Seaice Interface Variables === #include "ATM2D_VARS.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C dtatm, dtocn, dtcouple - Timesteps from couple.nml (hours) C myThid - Thread no. that called this routine. INTEGER dtatm, dtocn, dtcouple INTEGER myThid C LOCAL VARIABLES: INTEGER i,j,jj INTEGER ib, ibj1, ibj2 ! runoff band loop counters INTEGER j_atm, mn INTEGER dUnit _RL end1, end2, enda1, enda2, enda3 !used to compute grid conv areas _RL totrun_b(sNy) ! total file "runoff" in runoff bands _RL a1,a2 _RS atm_dyG(jm0) ! southern point/(boundary) of atmos grid DATA atm_dyG/2.0,44*4.0,2.0/ ! grid spacing for atmosphere dtatmo = dtatm * 3600. dtocno = dtocn * 3600. dtcouplo= dtcouple * 3600. C override data.ice seaice time step parms C these will need to change if coupling procedure changed thSice_deltaT = dtcouplo thsIce_dtTemp = dtatmo ocean_deltaT = dtcouplo CJRS This next check - only kill it if not MPI? IF (dtocno.NE.dTtracerLev(1)) THEN PRINT *,'Ocean tracer timestep differs between coupler ' PRINT *,'and the ocean data file' STOP ENDIF c Assuming the atmospheric grid array not passed, do this: atm_yG(1)=-90.0 DO j_atm=2,jm0 atm_yG(j_atm)=atm_yG(j_atm-1)+atm_dyG(j_atm-1) atm_yC(j_atm-1)=(atm_yG(j_atm-1)+atm_yG(j_atm))/2.0 ENDDO atm_yC(jm0)=atm_yG(jm0)+atm_dyG(jm0)/2.0 c end atmos grid initialization atm_oc_ind(1)=2 atm_oc_wgt(1)=1. _d 0 atm_oc_frac1(1)= (sin(yG(1,2,1,1)*deg2rad) - & sin(yG(1,1,1,1)*deg2rad))/ & (sin(atm_yG(3)*deg2rad)-sin(atm_yG(1)*deg2rad)) atm_oc_frac2(1)= 0. _d 0 ! assumes ocean(1) fits in atm(1) atm_oc_ind(sNy)=jm0-1 atm_oc_wgt(sNy)=1. _d 0 atm_oc_frac1(sNy)= (sin((yG(1,sNy,1,1) + & dyG(1,sNy,1,1)/6.37D6/deg2rad)*deg2rad)- & sin(yG(1,sNy,1,1)*deg2rad))/ & (sin((atm_yG(jm0)+atm_dyG(jm0))*deg2rad)- & sin(atm_yG(jm0-1)*deg2rad)) atm_oc_frac2(sNy)= 0. _d 0 ! assumes ocean(1) fits in atm(1) endwgt1 = sin(atm_yG(2)*deg2rad) !hard-coded that the atmos endwgt2 = sin(atm_yG(3)*deg2rad) - endwgt1 !grid is same in NH and SH endwgt1 = endwgt1 + 1. _d 0 !and goes 90S to 90N rsumwgt = 1. _d 0/(endwgt1 + endwgt2) atm_yG(2)=atm_yG(1) ! grid now combined atm end points atm_yG(jm0)=90. _d 0 DO j=2, sNy-1 DO jj=2,jm0-1 IF ((yG(1,j,1,1).GE.atm_yG(jj)).AND. & (yG(1,j,1,1).LT.atm_yG(jj+1))) j_atm=jj ENDDO atm_oc_ind(j)=j_atm end1= sin(yG(1,j,1,1) *deg2rad) end2= sin(yG(1,j+1,1,1) *deg2rad) enda1 = sin(atm_yG(j_atm) *deg2rad) enda2 = sin(atm_yG(j_atm+1) *deg2rad) IF ( yG(1,j+1,1,1) .GT. atm_yG(j_atm+1) ) THEN enda3 = sin(atm_yG(j_atm+2) *deg2rad) atm_oc_wgt(j)=(enda2-end1)/ (end2-end1) atm_oc_frac1(j)= (enda2-end1) / (enda2 - enda1) atm_oc_frac2(j)= (end2 - enda2) / (enda3 - enda2) ELSE atm_oc_wgt(j)=1. _d 0 atm_oc_frac1(j)= (end2-end1)/ (enda2-enda1) atm_oc_frac2(j)=0. _d 0 ENDIF ENDDO C compute tauv interpolation points tauv_jpt(1) = 2 ! south pole point; s/b land tauv_jwght(1) = 1. _d 0 DO j=2, sNy DO jj=1,jm0-1 IF (( yG(1,j,1,1) .GE. atm_yC(jj)).AND. & ( yG(1,j,1,1) .LT. atm_yC(jj+1))) j_atm=jj ENDDO tauv_jpt(j)= j_atm tauv_jwght(j)= 1. _d 0 - (yG(1,j,1,1) - atm_yC(j_atm)) / & (atm_yC(j_atm+1) - atm_yC(j_atm)) ENDDO C DO j=1,sNy C print *, 'j, tauv_jpt:', j,tauv_jpt(j),tauv_jwght(j) C ENDDO c c find land fraction c DO j_atm=1,jm0 cflan(j_atm)=0. _d 0 ocnArea(j_atm)=0. _d 0 ENDDO DO j=1,sNy DO i=1,sNx IF (maskC(i,j,1,1,1).EQ.1.) THEN ocnArea(atm_oc_ind(j))=ocnArea(atm_oc_ind(j)) + & rA(i,j,1,1)*atm_oc_wgt(j) IF (atm_oc_wgt(j).LT.1.d0) THEN ocnArea(atm_oc_ind(j)+1)=ocnArea(atm_oc_ind(j)+1) + & rA(i,j,1,1)*(1.d0-atm_oc_wgt(j)) ENDIF ENDIF ENDDO ENDDO DO j_atm=3,jm0-2 cflan(j_atm)=1. _d 0 - ocnArea(j_atm) / & (2. _d 0 * PI * 6.37 _d 6 * 6.37 _d 6 * & (sin(atm_yG(j_atm+1)*deg2rad) - sin(atm_yG(j_atm)*deg2rad))) if (cflan(j_atm).LT.1. _d -14) cflan(j_atm)=0. _d 0 ENDDO C deal with the combined atmos grid end cells... cflan(2)= 1. _d 0 - ocnArea(2) / & (2. _d 0*PI*6.37 _d 6*6.37 _d 6* & (sin(atm_yG(3)*deg2rad)+1. _d 0)) IF (cflan(2).LT.1. _d -14) cflan(2)=0. _d 0 cflan(1)=cflan(2) cflan(jm0-1)= 1.d0 - ocnArea(jm0-1) / & (2. _d 0*PI*6.37 _d 6*6.37 _d 6* & (1. _d 0-sin(atm_yG(jm0-1)*deg2rad))) IF (cflan(jm0-1).LT.1. _d -14) cflan(jm0-1)=0. _d 0 cflan(jm0)=cflan(jm0-1) PRINT *,'Land fractions on atmospheric grid: ' PRINT *, cflan PRINT *,'Lookup grid index, weights:' PRINT *, atm_oc_ind,atm_oc_wgt C PRINT *,'Lookup fraction 1 of atmos grid:' C PRINT *, atm_oc_frac1 C PRINT *,'Lookup fraction 2 of atmos grid:' C PRINT *, atm_oc_frac2 c c read in mean 1D atmos wind files -- store in memory c DO j_atm=1,jm0 DO mn=1,nForcingPer atau(j_atm,mn)=0. _d 0 atav(j_atm,mn)=0. _d 0 awind(j_atm,mn)=0. _d 0 ENDDO ENDDO CALL MDSFINDUNIT( dUnit, myThid ) IF ( atmosTauuFile .NE. ' ' ) THEN OPEN(dUnit, FILE=atmosTauuFile,STATUS='old', & ACCESS='direct', RECL=8*jm0*nForcingPer, & FORM='unformatted') READ(dUnit,REC=1), atau CLOSE(dUnit) ENDIF IF ( atmosTauvFile .NE. ' ' ) THEN OPEN(dUnit, FILE=atmosTauvFile, STATUS='old', & ACCESS='direct', RECL=8*jm0*nForcingPer, & FORM='unformatted') READ(dUnit, REC=1), atav CLOSE(dUnit) ENDIF IF ( atmosWindFile .NE. ' ' ) THEN OPEN(dUnit, FILE=atmosWindFile, STATUS='old', & ACCESS='direct', RECL=8*jm0*nForcingPer, & FORM='unformatted') READ(dUnit, REC=1), awind CLOSE(dUnit) ENDIF C The polar data point values for winds are effectively N/A given the C pole issue... although they are read in here, they are never used in C any calculations, as the polar ocean points access the data at atmos C 2 and jm0-1 points. c read in runoff data c to put runoff into specific grid cells c IF ( runoffMapFile .NE. ' ' ) THEN CALL READ_FLD_XY_RL( runoffMapFile, ' ', & runoffVal, 0, myThid ) ELSE DO j=1,sNy DO i=1,sNx if ( (maskC(i,j,1,1,1).EQ.1.) .AND. & ( (maskC(i-1,j,1,1,1).EQ.0.).OR. & (maskC(i+1,j,1,1,1).EQ.0.).OR. & (maskC(i,j-1,1,1,1).EQ.0.).OR. & (maskC(i,j+1,1,1,1).EQ.0.).OR. & (maskC(i+1,j+1,1,1,1).EQ.0.).OR. & (maskC(i-1,j-1,1,1,1).EQ.0.).OR. & (maskC(i+1,j-1,1,1,1).EQ.0.).OR. & (maskC(i-1,j+1,1,1,1).EQ.0.) ) ) THEN runoffVal(i,j)=1. _d 0 ENDIF ENDDO ENDDO ENDIF DO ib=1,numBands ibj1=1 if (ib.GT.1) ibj1= rband(ib-1)+1 ibj2=sNy if (ib.LT.numBands) ibj2= rband(ib) totrun_b(ib)=0.d0 DO j=ibj1,ibj2 DO i=1,sNx totrun_b(ib)=totrun_b(ib)+runoffVal(i,j)*maskC(i,j,1,1,1) ENDDO ENDDO DO j=ibj1,ibj2 runIndex(j)= ib ! for lookup of rband as fn. of latitude DO i=1,sNx runoffVal(i,j)=runoffVal(i,j)*maskC(i,j,1,1,1)/totrun_b(ib) ENDDO ENDDO ENDDO CALL INIT_SUMVARS(myThid) C Initialize 1D diagnostic variables DO j_atm=1,jm0 DO mn=1,nForcingPer sum_tauu_ta(j_atm,mn)= 0. _d 0 sum_tauv_ta(j_atm,mn)= 0. _d 0 sum_wsocean_ta(j_atm,mn)= 0. _d 0 sum_ps4ocean_ta(j_atm,mn)= 0. _d 0 ENDDO ENDDO C Initialize 2D diagnostic variables DO i=1-OLx,sNx+OLx DO j=1-OLy,sNy+OLy DO mn=1,nForcingPer qnet_atm_ta(i,j,mn)= 0. _d 0 evap_atm_ta(i,j,mn)= 0. _d 0 precip_atm_ta(i,j,mn)= 0. _d 0 runoff_atm_ta(i,j,mn)= 0. _d 0 sum_qrel_ta(i,j,mn)= 0. _d 0 sum_frel_ta(i,j,mn)= 0. _d 0 sum_iceMask_ta(i,j,mn)= 0. _d 0 sum_iceHeight_ta(i,j,mn)= 0. _d 0 sum_iceTime_ta(i,j,mn)= 0. _d 0 sum_oceMxLT_ta(i,j,mn)= 0. _d 0 sum_oceMxLS_ta(i,j,mn)= 0. _d 0 ENDDO qnet_atm(i,j)= 0. _d 0 evap_atm(i,j)= 0. _d 0 precip_atm(i,j)= 0. _d 0 runoff_atm(i,j)= 0. _d 0 sum_qrel(i,j)= 0. _d 0 sum_frel(i,j)= 0. _d 0 sum_iceMask(i,j)= 0. _d 0 sum_iceHeight(i,j)= 0. _d 0 sum_iceTime(i,j)= 0. _d 0 sum_oceMxLT(i,j)= 0. _d 0 sum_oceMxLS(i,j)= 0. _d 0 ENDDO ENDDO C Initialize year-end diags and max/min seaice variables SHice_min = 1. _d 18 NHice_min = 1. _d 18 SHice_max = 0. _d 0 NHice_max = 0. _d 0 sst_tave= 0. _d 0 sss_tave= 0. _d 0 HF2ocn_tave= 0. _d 0 FW2ocn_tave= 0. _d 0 CO2flx_tave= 0. _d 0 OPEN(25,FILE='resocean.dat',STATUS='replace') CLOSE(25) C Initialize following for safety and/or cold start DO i=1-OLx,sNx+OLx DO j=1-OLy,sNy+OLy pass_runoff(i,j)= 0. _d 0 pass_qnet(i,j)= 0. _d 0 pass_evap(i,j)= 0. _d 0 pass_precip(i,j)= 0. _d 0 pass_fu(i,j)= 0. _d 0 pass_fv(i,j)= 0. _d 0 pass_wspeed(i,j)= 0. _d 0 pass_solarnet(i,j)= 0. _d 0 pass_slp(i,j)= 0. _d 0 pass_siceLoad(i,j)= 0. _d 0 pass_pCO2(i,j)= 0. _d 0 pass_prcAtm(i,j) = 0. _d 0 snowPrc (i,j) = 0. _d 0 sFluxFromIce(i,j)= 0. _d 0 ENDDO ENDDO C Initialize following (if ocn carbon not passed) DO i=1,sNx DO j=1,sNy oFluxCO2(i,j) = 0. _d 0 ENDDO ENDDO RETURN END