| 1 |
C $Header: /u/gcmpack/MITgcm/pkg/atm2d/init_atm2d.F,v 1.6 2007/11/19 22:57:31 jscott Exp $ |
| 2 |
C $Name: $ |
| 3 |
|
| 4 |
#include "ctrparam.h" |
| 5 |
#include "ATM2D_OPTIONS.h" |
| 6 |
|
| 7 |
C !INTERFACE: |
| 8 |
SUBROUTINE INIT_ATM2D(dtatm, dtocn, dtcouple, myThid ) |
| 9 |
C *==========================================================* |
| 10 |
C | INIT_1DTO2D | |
| 11 |
C | This initialization routine should be run after the | |
| 12 |
c | the ocean grid/pickup have been read in. | |
| 13 |
c | | |
| 14 |
c | Note: grid variable indices bi,bj are hard-coded 1,1 | |
| 15 |
c | This should work if coupler or atmos/coupler on one | |
| 16 |
c | machine. | |
| 17 |
c | | |
| 18 |
C *==========================================================* |
| 19 |
c |
| 20 |
IMPLICIT NONE |
| 21 |
|
| 22 |
C === Global Atmosphere Variables === |
| 23 |
#include "ATMSIZE.h" |
| 24 |
#include "AGRID.h" |
| 25 |
|
| 26 |
C === Global Ocean Variables === |
| 27 |
#include "SIZE.h" |
| 28 |
#include "EEPARAMS.h" |
| 29 |
#include "PARAMS.h" |
| 30 |
#include "GRID.h" |
| 31 |
|
| 32 |
C === Global SeaIce Parameters === |
| 33 |
#include "THSICE_PARAMS.h" |
| 34 |
|
| 35 |
C === Atmos/Ocean/Seaice Interface Variables === |
| 36 |
#include "ATM2D_VARS.h" |
| 37 |
|
| 38 |
|
| 39 |
C !INPUT/OUTPUT PARAMETERS: |
| 40 |
C === Routine arguments === |
| 41 |
C dtatm, dtocn, dtcouple - Timesteps from couple.nml (hours) |
| 42 |
C myThid - Thread no. that called this routine. |
| 43 |
INTEGER dtatm, dtocn, dtcouple |
| 44 |
INTEGER myThid |
| 45 |
|
| 46 |
C LOCAL VARIABLES: |
| 47 |
INTEGER i,j,jj |
| 48 |
INTEGER ib, ibj1, ibj2 ! runoff band loop counters |
| 49 |
INTEGER j_atm, mn |
| 50 |
INTEGER dUnit |
| 51 |
_RL end1, end2, enda1, enda2, enda3 !used to compute grid conv areas |
| 52 |
_RL totrun_b(sNy) ! total file "runoff" in runoff bands |
| 53 |
_RL a1,a2 |
| 54 |
_RS atm_dyG(jm0) ! southern point/(boundary) of atmos grid |
| 55 |
DATA atm_dyG/2.0,44*4.0,2.0/ ! grid spacing for atmosphere |
| 56 |
|
| 57 |
dtatmo = dtatm * 3600. |
| 58 |
dtocno = dtocn * 3600. |
| 59 |
dtcouplo= dtcouple * 3600. |
| 60 |
|
| 61 |
C override data.ice seaice time step parms |
| 62 |
C these will need to change if coupling procedure changed |
| 63 |
thSice_deltaT = dtcouplo |
| 64 |
thsIce_dtTemp = dtatmo |
| 65 |
ocean_deltaT = dtcouplo |
| 66 |
|
| 67 |
CJRS This next check - only kill it if not MPI? |
| 68 |
IF (dtocno.NE.dTtracerLev(1)) THEN |
| 69 |
PRINT *,'Ocean tracer timestep differs between coupler ' |
| 70 |
PRINT *,'and the ocean data file' |
| 71 |
STOP |
| 72 |
ENDIF |
| 73 |
|
| 74 |
c Assuming the atmospheric grid array not passed, do this: |
| 75 |
atm_yG(1)=-90.0 |
| 76 |
DO j_atm=2,jm0 |
| 77 |
atm_yG(j_atm)=atm_yG(j_atm-1)+atm_dyG(j_atm-1) |
| 78 |
atm_yC(j_atm-1)=(atm_yG(j_atm-1)+atm_yG(j_atm))/2.0 |
| 79 |
ENDDO |
| 80 |
atm_yC(jm0)=atm_yG(jm0)+atm_dyG(jm0)/2.0 |
| 81 |
|
| 82 |
c end atmos grid initialization |
| 83 |
|
| 84 |
atm_oc_ind(1)=2 |
| 85 |
atm_oc_wgt(1)=1. _d 0 |
| 86 |
atm_oc_frac1(1)= (sin(yG(1,2,1,1)*deg2rad) - |
| 87 |
& sin(yG(1,1,1,1)*deg2rad))/ |
| 88 |
& (sin(atm_yG(3)*deg2rad)-sin(atm_yG(1)*deg2rad)) |
| 89 |
atm_oc_frac2(1)= 0. _d 0 ! assumes ocean(1) fits in atm(1) |
| 90 |
atm_oc_ind(sNy)=jm0-1 |
| 91 |
atm_oc_wgt(sNy)=1. _d 0 |
| 92 |
atm_oc_frac1(sNy)= (sin((yG(1,sNy,1,1) + |
| 93 |
& dyG(1,sNy,1,1)/6.37D6/deg2rad)*deg2rad)- |
| 94 |
& sin(yG(1,sNy,1,1)*deg2rad))/ |
| 95 |
& (sin((atm_yG(jm0)+atm_dyG(jm0))*deg2rad)- |
| 96 |
& sin(atm_yG(jm0-1)*deg2rad)) |
| 97 |
atm_oc_frac2(sNy)= 0. _d 0 ! assumes ocean(1) fits in atm(1) |
| 98 |
|
| 99 |
endwgt1 = sin(atm_yG(2)*deg2rad) !hard-coded that the atmos |
| 100 |
endwgt2 = sin(atm_yG(3)*deg2rad) - endwgt1 !grid is same in NH and SH |
| 101 |
endwgt1 = endwgt1 + 1. _d 0 !and goes 90S to 90N |
| 102 |
rsumwgt = 1. _d 0/(endwgt1 + endwgt2) |
| 103 |
|
| 104 |
atm_yG(2)=atm_yG(1) ! grid now combined atm end points |
| 105 |
atm_yG(jm0)=90. _d 0 |
| 106 |
|
| 107 |
DO j=2, sNy-1 |
| 108 |
|
| 109 |
DO jj=2,jm0-1 |
| 110 |
IF ((yG(1,j,1,1).GE.atm_yG(jj)).AND. |
| 111 |
& (yG(1,j,1,1).LT.atm_yG(jj+1))) j_atm=jj |
| 112 |
ENDDO |
| 113 |
|
| 114 |
atm_oc_ind(j)=j_atm |
| 115 |
end1= sin(yG(1,j,1,1) *deg2rad) |
| 116 |
end2= sin(yG(1,j+1,1,1) *deg2rad) |
| 117 |
enda1 = sin(atm_yG(j_atm) *deg2rad) |
| 118 |
enda2 = sin(atm_yG(j_atm+1) *deg2rad) |
| 119 |
IF ( yG(1,j+1,1,1) .GT. atm_yG(j_atm+1) ) THEN |
| 120 |
enda3 = sin(atm_yG(j_atm+2) *deg2rad) |
| 121 |
atm_oc_wgt(j)=(enda2-end1)/ (end2-end1) |
| 122 |
atm_oc_frac1(j)= (enda2-end1) / (enda2 - enda1) |
| 123 |
atm_oc_frac2(j)= (end2 - enda2) / (enda3 - enda2) |
| 124 |
ELSE |
| 125 |
atm_oc_wgt(j)=1. _d 0 |
| 126 |
atm_oc_frac1(j)= (end2-end1)/ (enda2-enda1) |
| 127 |
atm_oc_frac2(j)=0. _d 0 |
| 128 |
ENDIF |
| 129 |
ENDDO |
| 130 |
|
| 131 |
C compute tauv interpolation points |
| 132 |
tauv_jpt(1) = 2 ! south pole point; s/b land |
| 133 |
tauv_jwght(1) = 1. _d 0 |
| 134 |
DO j=2, sNy |
| 135 |
DO jj=1,jm0-1 |
| 136 |
IF (( yG(1,j,1,1) .GE. atm_yC(jj)).AND. |
| 137 |
& ( yG(1,j,1,1) .LT. atm_yC(jj+1))) j_atm=jj |
| 138 |
ENDDO |
| 139 |
tauv_jpt(j)= j_atm |
| 140 |
tauv_jwght(j)= 1. _d 0 - (yG(1,j,1,1) - atm_yC(j_atm)) / |
| 141 |
& (atm_yC(j_atm+1) - atm_yC(j_atm)) |
| 142 |
ENDDO |
| 143 |
|
| 144 |
C DO j=1,sNy |
| 145 |
C print *, 'j, tauv_jpt:', j,tauv_jpt(j),tauv_jwght(j) |
| 146 |
C ENDDO |
| 147 |
|
| 148 |
c |
| 149 |
c find land fraction |
| 150 |
c |
| 151 |
DO j_atm=1,jm0 |
| 152 |
cflan(j_atm)=0. _d 0 |
| 153 |
ocnArea(j_atm)=0. _d 0 |
| 154 |
ENDDO |
| 155 |
|
| 156 |
DO j=1,sNy |
| 157 |
DO i=1,sNx |
| 158 |
IF (maskC(i,j,1,1,1).EQ.1.) THEN |
| 159 |
ocnArea(atm_oc_ind(j))=ocnArea(atm_oc_ind(j)) + |
| 160 |
& rA(i,j,1,1)*atm_oc_wgt(j) |
| 161 |
IF (atm_oc_wgt(j).LT.1.d0) THEN |
| 162 |
ocnArea(atm_oc_ind(j)+1)=ocnArea(atm_oc_ind(j)+1) + |
| 163 |
& rA(i,j,1,1)*(1.d0-atm_oc_wgt(j)) |
| 164 |
ENDIF |
| 165 |
ENDIF |
| 166 |
ENDDO |
| 167 |
ENDDO |
| 168 |
|
| 169 |
DO j_atm=3,jm0-2 |
| 170 |
cflan(j_atm)=1. _d 0 - ocnArea(j_atm) / |
| 171 |
& (2. _d 0 * PI * 6.37 _d 6 * 6.37 _d 6 * |
| 172 |
& (sin(atm_yG(j_atm+1)*deg2rad) - sin(atm_yG(j_atm)*deg2rad))) |
| 173 |
if (cflan(j_atm).LT.1. _d -14) cflan(j_atm)=0. _d 0 |
| 174 |
ENDDO |
| 175 |
|
| 176 |
C deal with the combined atmos grid end cells... |
| 177 |
cflan(2)= 1. _d 0 - ocnArea(2) / |
| 178 |
& (2. _d 0*PI*6.37 _d 6*6.37 _d 6* |
| 179 |
& (sin(atm_yG(3)*deg2rad)+1. _d 0)) |
| 180 |
IF (cflan(2).LT.1. _d -14) cflan(2)=0. _d 0 |
| 181 |
cflan(1)=cflan(2) |
| 182 |
cflan(jm0-1)= 1.d0 - ocnArea(jm0-1) / |
| 183 |
& (2. _d 0*PI*6.37 _d 6*6.37 _d 6* |
| 184 |
& (1. _d 0-sin(atm_yG(jm0-1)*deg2rad))) |
| 185 |
IF (cflan(jm0-1).LT.1. _d -14) cflan(jm0-1)=0. _d 0 |
| 186 |
cflan(jm0)=cflan(jm0-1) |
| 187 |
|
| 188 |
PRINT *,'Land fractions on atmospheric grid: ' |
| 189 |
PRINT *, cflan |
| 190 |
PRINT *,'Lookup grid index, weights:' |
| 191 |
PRINT *, atm_oc_ind,atm_oc_wgt |
| 192 |
C PRINT *,'Lookup fraction 1 of atmos grid:' |
| 193 |
C PRINT *, atm_oc_frac1 |
| 194 |
C PRINT *,'Lookup fraction 2 of atmos grid:' |
| 195 |
C PRINT *, atm_oc_frac2 |
| 196 |
|
| 197 |
c |
| 198 |
c read in mean 1D atmos wind files -- store in memory |
| 199 |
c |
| 200 |
DO j_atm=1,jm0 |
| 201 |
DO mn=1,nForcingPer |
| 202 |
atau(j_atm,mn)=0. _d 0 |
| 203 |
atav(j_atm,mn)=0. _d 0 |
| 204 |
awind(j_atm,mn)=0. _d 0 |
| 205 |
ENDDO |
| 206 |
ENDDO |
| 207 |
|
| 208 |
CALL MDSFINDUNIT( dUnit, myThid ) |
| 209 |
|
| 210 |
IF ( atmosTauuFile .NE. ' ' ) THEN |
| 211 |
OPEN(dUnit, FILE=atmosTauuFile,STATUS='old', |
| 212 |
& ACCESS='direct', RECL=8*jm0*nForcingPer, |
| 213 |
& FORM='unformatted') |
| 214 |
READ(dUnit,REC=1), atau |
| 215 |
CLOSE(dUnit) |
| 216 |
ENDIF |
| 217 |
|
| 218 |
IF ( atmosTauvFile .NE. ' ' ) THEN |
| 219 |
OPEN(dUnit, FILE=atmosTauvFile, STATUS='old', |
| 220 |
& ACCESS='direct', RECL=8*jm0*nForcingPer, |
| 221 |
& FORM='unformatted') |
| 222 |
READ(dUnit, REC=1), atav |
| 223 |
CLOSE(dUnit) |
| 224 |
ENDIF |
| 225 |
|
| 226 |
IF ( atmosWindFile .NE. ' ' ) THEN |
| 227 |
OPEN(dUnit, FILE=atmosWindFile, STATUS='old', |
| 228 |
& ACCESS='direct', RECL=8*jm0*nForcingPer, |
| 229 |
& FORM='unformatted') |
| 230 |
READ(dUnit, REC=1), awind |
| 231 |
CLOSE(dUnit) |
| 232 |
ENDIF |
| 233 |
|
| 234 |
C The polar data point values for winds are effectively N/A given the |
| 235 |
C pole issue... although they are read in here, they are never used in |
| 236 |
C any calculations, as the polar ocean points access the data at atmos |
| 237 |
C 2 and jm0-1 points. |
| 238 |
|
| 239 |
|
| 240 |
c read in runoff data |
| 241 |
c to put runoff into specific grid cells |
| 242 |
c |
| 243 |
IF ( runoffMapFile .NE. ' ' ) THEN |
| 244 |
CALL READ_FLD_XY_RL( runoffMapFile, ' ', |
| 245 |
& runoffVal, 0, myThid ) |
| 246 |
ELSE |
| 247 |
DO j=1,sNy |
| 248 |
DO i=1,sNx |
| 249 |
if ( (maskC(i,j,1,1,1).EQ.1.) .AND. |
| 250 |
& ( (maskC(i-1,j,1,1,1).EQ.0.).OR. |
| 251 |
& (maskC(i+1,j,1,1,1).EQ.0.).OR. |
| 252 |
& (maskC(i,j-1,1,1,1).EQ.0.).OR. |
| 253 |
& (maskC(i,j+1,1,1,1).EQ.0.).OR. |
| 254 |
& (maskC(i+1,j+1,1,1,1).EQ.0.).OR. |
| 255 |
& (maskC(i-1,j-1,1,1,1).EQ.0.).OR. |
| 256 |
& (maskC(i+1,j-1,1,1,1).EQ.0.).OR. |
| 257 |
& (maskC(i-1,j+1,1,1,1).EQ.0.) ) ) THEN |
| 258 |
runoffVal(i,j)=1. _d 0 |
| 259 |
ENDIF |
| 260 |
ENDDO |
| 261 |
ENDDO |
| 262 |
ENDIF |
| 263 |
|
| 264 |
DO ib=1,numBands |
| 265 |
ibj1=1 |
| 266 |
if (ib.GT.1) ibj1= rband(ib-1)+1 |
| 267 |
ibj2=sNy |
| 268 |
if (ib.LT.numBands) ibj2= rband(ib) |
| 269 |
totrun_b(ib)=0.d0 |
| 270 |
DO j=ibj1,ibj2 |
| 271 |
DO i=1,sNx |
| 272 |
totrun_b(ib)=totrun_b(ib)+runoffVal(i,j)*maskC(i,j,1,1,1) |
| 273 |
ENDDO |
| 274 |
ENDDO |
| 275 |
DO j=ibj1,ibj2 |
| 276 |
runIndex(j)= ib ! for lookup of rband as fn. of latitude |
| 277 |
DO i=1,sNx |
| 278 |
runoffVal(i,j)=runoffVal(i,j)*maskC(i,j,1,1,1)/totrun_b(ib) |
| 279 |
ENDDO |
| 280 |
ENDDO |
| 281 |
ENDDO |
| 282 |
|
| 283 |
CALL INIT_SUMVARS(myThid) |
| 284 |
|
| 285 |
C Initialize 1D diagnostic variables |
| 286 |
DO j_atm=1,jm0 |
| 287 |
DO mn=1,nForcingPer |
| 288 |
sum_tauu_ta(j_atm,mn)= 0. _d 0 |
| 289 |
sum_tauv_ta(j_atm,mn)= 0. _d 0 |
| 290 |
sum_wsocean_ta(j_atm,mn)= 0. _d 0 |
| 291 |
sum_ps4ocean_ta(j_atm,mn)= 0. _d 0 |
| 292 |
ENDDO |
| 293 |
ENDDO |
| 294 |
|
| 295 |
C Initialize 2D diagnostic variables |
| 296 |
DO i=1-OLx,sNx+OLx |
| 297 |
DO j=1-OLy,sNy+OLy |
| 298 |
DO mn=1,nForcingPer |
| 299 |
qnet_atm_ta(i,j,mn)= 0. _d 0 |
| 300 |
evap_atm_ta(i,j,mn)= 0. _d 0 |
| 301 |
precip_atm_ta(i,j,mn)= 0. _d 0 |
| 302 |
runoff_atm_ta(i,j,mn)= 0. _d 0 |
| 303 |
sum_qrel_ta(i,j,mn)= 0. _d 0 |
| 304 |
sum_frel_ta(i,j,mn)= 0. _d 0 |
| 305 |
sum_iceMask_ta(i,j,mn)= 0. _d 0 |
| 306 |
sum_iceHeight_ta(i,j,mn)= 0. _d 0 |
| 307 |
sum_iceTime_ta(i,j,mn)= 0. _d 0 |
| 308 |
sum_oceMxLT_ta(i,j,mn)= 0. _d 0 |
| 309 |
sum_oceMxLS_ta(i,j,mn)= 0. _d 0 |
| 310 |
ENDDO |
| 311 |
qnet_atm(i,j)= 0. _d 0 |
| 312 |
evap_atm(i,j)= 0. _d 0 |
| 313 |
precip_atm(i,j)= 0. _d 0 |
| 314 |
runoff_atm(i,j)= 0. _d 0 |
| 315 |
sum_qrel(i,j)= 0. _d 0 |
| 316 |
sum_frel(i,j)= 0. _d 0 |
| 317 |
sum_iceMask(i,j)= 0. _d 0 |
| 318 |
sum_iceHeight(i,j)= 0. _d 0 |
| 319 |
sum_iceTime(i,j)= 0. _d 0 |
| 320 |
sum_oceMxLT(i,j)= 0. _d 0 |
| 321 |
sum_oceMxLS(i,j)= 0. _d 0 |
| 322 |
ENDDO |
| 323 |
ENDDO |
| 324 |
|
| 325 |
C Initialize year-end diags and max/min seaice variables |
| 326 |
SHice_min = 1. _d 18 |
| 327 |
NHice_min = 1. _d 18 |
| 328 |
SHice_max = 0. _d 0 |
| 329 |
NHice_max = 0. _d 0 |
| 330 |
sst_tave= 0. _d 0 |
| 331 |
sss_tave= 0. _d 0 |
| 332 |
HF2ocn_tave= 0. _d 0 |
| 333 |
FW2ocn_tave= 0. _d 0 |
| 334 |
CO2flx_tave= 0. _d 0 |
| 335 |
OPEN(25,FILE='resocean.dat',STATUS='replace') |
| 336 |
CLOSE(25) |
| 337 |
|
| 338 |
C Initialize following for safety and/or cold start |
| 339 |
DO i=1-OLx,sNx+OLx |
| 340 |
DO j=1-OLy,sNy+OLy |
| 341 |
pass_runoff(i,j)= 0. _d 0 |
| 342 |
pass_qnet(i,j)= 0. _d 0 |
| 343 |
pass_evap(i,j)= 0. _d 0 |
| 344 |
pass_precip(i,j)= 0. _d 0 |
| 345 |
pass_fu(i,j)= 0. _d 0 |
| 346 |
pass_fv(i,j)= 0. _d 0 |
| 347 |
pass_wspeed(i,j)= 0. _d 0 |
| 348 |
pass_solarnet(i,j)= 0. _d 0 |
| 349 |
pass_slp(i,j)= 0. _d 0 |
| 350 |
pass_siceLoad(i,j)= 0. _d 0 |
| 351 |
pass_pCO2(i,j)= 0. _d 0 |
| 352 |
pass_prcAtm(i,j)= 0. _d 0 |
| 353 |
sFluxFromIce(i,j)= 0. _d 0 |
| 354 |
ENDDO |
| 355 |
ENDDO |
| 356 |
|
| 357 |
C Initialize following (if ocn carbon not passed) |
| 358 |
DO i=1,sNx |
| 359 |
DO j=1,sNy |
| 360 |
oFluxCO2(i,j) = 0. _d 0 |
| 361 |
ENDDO |
| 362 |
ENDDO |
| 363 |
|
| 364 |
RETURN |
| 365 |
END |
| 366 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph