pkg/atm2d/forward_step_atm2d.F

09a6f3668a Jeff*0001 #include "ctrparam.h"
                0002 #ifdef OCEAN_3D
                0003 #  include "ATM2D_OPTIONS.h"
                0004 #endif
                0005 C
                0006       SUBROUTINE FORWARD_STEP_ATM2D(iloop, myTime, myIter, myThid)
                0007 C     |==========================================================|
                0008 C     | Does time loop for one coupled period. The main loop     |
                0009 C     | this is the MITGCM main loop OR a separate driver for    |
                0010 C     | IGSM 2.2                                                 |
                0011 C     \==========================================================/
                0012       IMPLICIT NONE
                0013
                0014 #include "ATMSIZE.h"
                0015 #include "DRIVER.h"
                0016
                0017 #ifdef OCEAN_3D
                0018 #  include "SIZE.h"
                0019 #  include "EEPARAMS.h"
                0020 #  include "PARAMS.h"
                0021 #  include "ATM2D_VARS.h"
                0022 #endif
                0023
50a1736d54 Jean*0024 #ifdef NCEPWIND
260186e531 Jeff*0025       COMMON  /SEED/JSEED,IFRST,NEXTN
                0026       INTEGER JSEED,IFRST,NEXTN
                0027 #endif
                0028
09a6f3668a Jeff*0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     == Routine arguments ==
                0031 C     iloop  - loop counter for coupled period time steps (main time loop)
                0032 C     myIter - iteration counter for this thread (ocean times steps +nIter0)
                0033 C     myTime - time counter for this thread (ocean time, from starttTime)
                0034 C     myThid - thread number for this instance of the routine.
                0035       INTEGER iloop
                0036       REAL*8  myTime
                0037       INTEGER myIter
                0038       INTEGER myThid
                0039
                0040 C     === Local variables ===
                0041       INTEGER idyear ! year # of simulation, starting at year 1
                0042       INTEGER iyr    ! year # of simulation, starting from specified inyear
                0043       INTEGER inyr   ! hours into the current year, end of coupled period
                0044       INTEGER monid  ! current month of the year
                0045       INTEGER inday  ! hour of the day, end of the coupled period
                0046       INTEGER dayid  ! day of the current month
                0047       INTEGER j,mn,na,no   !loop counters
                0048       INTEGER jdofmhr(0:12)
                0049       DATA jdofmhr/0,744,1416,2160,2880,3624,4344,5088,
                0050      &                  5832,6552,7296,8016,8760/
                0051 C i.e. 0,31*24,59*24,90*24,120*24,151*24,181*24,
                0052 C      212*24,243*24,273*24,304*24,334*24,365*24
                0053 #ifdef CPL_TEM
                0054       INTEGER ndmonth(12)
                0055       DATA ndmonth/31,28,31,30,31,30,31,31,30,31,30,31/
                0056       REAL*4 totup, aduptt
084d5af5e3 Jeff*0057       REAL*8 tcumn
c6a279aff5 Jeff*0058 #endif
                0059 #if (defined CPL_TEM) || (defined CPL_OCEANCO2)
5581185843 Jeff*0060       REAL*8 nepan
50a1736d54 Jean*0061       REAL*8 ocuptp
09a6f3668a Jeff*0062 #endif
                0063 #ifdef OCEAN_3D
50a1736d54 Jean*0064       INTEGER iloop_ocn, i
                0065       _RL qPrcRn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0066 #  ifdef NCEPWIND
260186e531 Jeff*0067       REAL*8  RAND
                0068       CHARACTER *4 ncep_yr
                0069 #  endif
09a6f3668a Jeff*0070 #endif
50a1736d54 Jean*0071 #ifdef DATA4TEM
260186e531 Jeff*0072       CHARACTER *40 f4tem,f4clm,f24tem,f34tem
5581185843 Jeff*0073       CHARACTER *4  cfile
                0074       CHARACTER *8  f14tem,f14clm
260186e531 Jeff*0075       character *9 f124tem,f134tem
5581185843 Jeff*0076       f14tem='data4tem'
                0077       f14clm='data4clm'
260186e531 Jeff*0078       f124tem='data24tem'
                0079       f134tem='data34tem'
                0080       nfile=1
5581185843 Jeff*0081 #endif
09a6f3668a Jeff*0082
f345ec4dc7 Jeff*0083 C     print *,'*** Top of forwrdstep_atm',iloop,myTime,myIter
09a6f3668a Jeff*0084       idyear= int((iloop-1)*dtcouple/365.0/24.0) + 1
                0085       iyr= idyear + startYear -1
                0086       inyr = mod(iloop*dtcouple, 365*24)
8ad2868148 Jeff*0087       IF (inyr .EQ. 0) inyr=jdofmhr(12)
09a6f3668a Jeff*0088       DO mn=1,12
                0089         IF ((inyr.GT.jdofmhr(mn-1)).AND.(inyr.LE.jdofmhr(mn))) monid=mn
                0090       ENDDO
                0091       inday= mod(iloop*dtcouple, 24)
                0092       dayid= int((inyr-dtcouple-jdofmhr(monid-1))/24.0) +1
f345ec4dc7 Jeff*0093 C     print *,'*** idyear,iyr,inyr,monid,inday,dayid',
                0094 C    &             idyear,iyr,inyr,monid,inday,dayid
09a6f3668a Jeff*0095
                0096       IF (inyr.EQ.dtcouple) THEN !do this block at start of new year
                0097         PRINT *,'*** Starting a new year'
50a1736d54 Jean*0098 #ifdef NCEPWIND
260186e531 Jeff*0099         WRITE(ncep_yr,'(I4)') (NINT(RAND()*60.0+0.5) + 1947)
                0100         PRINT *,'Using NCEP wind variations from year: ',ncep_yr
50a1736d54 Jean*0101          OPEN(6007,
260186e531 Jeff*0102      &     FILE='ncep_taux_variations_'//ncep_yr//'.bin',STATUS='old',
                0103      &     ACCESS='direct', RECL=4*sNx*sNy,
                0104      &     FORM='unformatted')
50a1736d54 Jean*0105          OPEN(6008,
260186e531 Jeff*0106      &     FILE='ncep_tauy_variations_'//ncep_yr//'.bin',STATUS='old',
                0107      &     ACCESS='direct', RECL=4*sNx*sNy,
                0108      &     FORM='unformatted')
50a1736d54 Jean*0109          OPEN(6009,
260186e531 Jeff*0110      &     FILE='ncep_speed_variations_'//ncep_yr//'.bin',STATUS='old',
                0111      &     ACCESS='direct', RECL=4*sNx*sNy,
                0112      &     FORM='unformatted')
                0113          ncep_counter = 1
                0114 #endif
09a6f3668a Jeff*0115 #ifdef DATA4TEM
                0116         IF (nfile.gt.1)THEN
                0117           CLOSE(935)
                0118           CLOSE(937)
260186e531 Jeff*0119           CLOSE(938)
                0120           CLOSE(939)
09a6f3668a Jeff*0121         ENDIF
                0122         IF(iyr.gt.1000) THEN
                0123            nfile=iyr
                0124         ELSE
                0125           nfile=1000+iyr
                0126         ENDIF
                0127         WRITE (cfile,'i4'),nfile
                0128         f4tem=f14tem//cfile
                0129         f4clm=f14clm//cfile
260186e531 Jeff*0130         f24tem=f124tem//cfile
                0131         f34tem=f134tem//cfile
09a6f3668a Jeff*0132         OPEN(935,file=f4clm,form='unformatted',status='new')
                0133         OPEN(937,file=f4tem,form='unformatted',status='new')
260186e531 Jeff*0134         OPEN(938,file=f24tem,form='unformatted',status='new')
                0135         OPEN(939,file=f34tem,form='unformatted',status='new')
09a6f3668a Jeff*0136         nfile=nfile+1
                0137 #endif
                0138 #ifdef CPL_TEM
                0139         nepan=0.0
                0140         ch4ann=0.0
                0141         n2oann=0.0
                0142         xco2ann=0.0
                0143 #endif
                0144 #ifdef CPL_OCEANCO2
                0145         temuptann=0.
                0146         DO j=1,jm0
                0147           co24ocnan(j)=0.0
                0148         ENDDO
d09f90e9d3 Jeff*0149 #  ifdef ML_2D
49512b0964 Jeff*0150           call kvcarbon(iyr)
d09f90e9d3 Jeff*0151 #  endif
09a6f3668a Jeff*0152 #endif
                0153 #ifdef CPL_TEM
                0154         DO j=1,jm0
                0155           antemnep(j)=0.
                0156         ENDDO
                0157 #  ifndef CPL_CHEM
                0158         CALL robso3(iyr)
                0159 #  endif
                0160 C For land use
                0161         CALL updatelcluc(idyear)
                0162 #endif
7090ec812e Jeff*0163 #ifdef CPL_CHEM
                0164         print *,' Before eppaemission'
                0165         CALL eppaemission (iyr)
                0166 #endif
09a6f3668a Jeff*0167       ENDIF   !end block done at year-start
                0168
                0169       IF (inyr.EQ.jdofmhr(monid-1)+dtcouple) THEN !do this block month start
                0170        PRINT *,'***Starting a new month'
                0171 #ifdef CPL_TEM
                0172         CALL zclimate2tem
                0173 #endif
                0174 #ifdef CPL_OCEANCO2
                0175         ocumn=0.0
                0176         DO j=1,jm0
                0177           fluxco2mn(j)=0.0
                0178         ENDDO
                0179 #endif
7ad882c4ff Jeff*0180 #ifdef OCEAN_3D
                0181         new_mon= .TRUE.
                0182 #endif
09a6f3668a Jeff*0183       ENDIF  !end block at start of the month
                0184 C
                0185 C------------------- Top of Coupled Period Loop --------------------------
                0186 C
                0187
                0188 #ifdef OCEAN_3D
260186e531 Jeff*0189 #  ifdef NCEPWIND
                0190 C     Read in the next timestep of ncep wind stress variations
                0191 C      PRINT *,'*** Read in next NCEP record'
                0192       READ(6007, REC=ncep_counter), fu_ncep
                0193       READ(6008, REC=ncep_counter), fv_ncep
                0194       READ(6009, REC=ncep_counter), fs_ncep
                0195       ncep_counter=ncep_counter+1
                0196 #  endif
9274434acc Jean*0197 #  ifdef ATM2D_MPI_ON
09a6f3668a Jeff*0198       CALL CPL_RECV_OCN_FIELDS
                0199 #  endif
                0200       CALL GET_OCNVARS( myTime, myIter, myThid)
                0201       IF ( (iloop.NE.1).OR. (iloop.EQ.1.AND.
                0202      &      (startTime.NE.baseTime .OR. nIter0.NE.0)) ) THEN
                0203 C       don't run the ice growth/melt on step 1 if "cold" start
50a1736d54 Jean*0204         DO j = 1-OLy, sNy+OLy
                0205          DO i = 1-OLx, sNx+OLx
                0206           qPrcRn(i,j) = 0.
                0207          ENDDO
                0208         ENDDO
7b7831d041 Jean*0209         CALL THSICE_STEP_FWD( 1, 1, 1, sNx, 1, sNy,
                0210      I                        pass_prcAtm, snowPrc, qPrcRn,
                0211      &                        myTime, myIter, myThid )
09a6f3668a Jeff*0212         CALL THSICE_AVE( 1,1, myTime, myIter, myThid )
                0213       ENDIF
                0214       CALL CALC_ZONAL_MEANS(.TRUE.,myThid)
                0215       CALL PUT_OCNVARS(myTime,myIter,myThid)
8e101fde6e Jeff*0216       CALL SUM_YR_END_DIAGS(myTime,myIter,myThid)
09a6f3668a Jeff*0217 #  ifdef ATM2D_MPI_ON
                0218       CALL CPL_SEND_OCN_FIELDS
                0219 #  endif
9274434acc Jean*0220 #endif
                0221
b6a1ae81d9 Jeff*0222 C      PRINT *,'Top of ncall_atm Loop'
09a6f3668a Jeff*0223       DO na=1,ncall_atm    !loop for atmos forward time steps
                0224         CALL atmosphere(dtatm,monid)
                0225 #ifdef OCEAN_3D
                0226         CALL ATM2OCN_MAIN(iloop, na, monid, myIter, myThid)
                0227         CALL SUM_OCN_FLUXES(myThid)
aea46876b4 Jeff*0228         CALL PASS_THSICE_FLUXES(myThid)
b6a1ae81d9 Jeff*0229         CALL THSICE_IMPL_TEMP(netSW, sFlx, dTsurf, 1,1,
09a6f3668a Jeff*0230      &                        myTime, myIter, myThid)
aea46876b4 Jeff*0231         CALL SUM_THSICE_OUT(myThid)
09a6f3668a Jeff*0232         CALL CALC_ZONAL_MEANS(.FALSE.,myThid) !just mean Tsrf recalculated
                0233 #endif
                0234       ENDDO  ! ncall_atm loop
                0235
b6a1ae81d9 Jeff*0236 C      PRINT *,'Top of ncall_ocean Loop'
09a6f3668a Jeff*0237       DO no=1,ncall_ocean   !loop for each ocean forward step
                0238
9274434acc Jean*0239 #ifdef OCEAN_3D
09a6f3668a Jeff*0240         iloop_ocn = nint((iloop-1)*dtcouplo/deltaTClock) + no
                0241 #  ifndef ATM2D_MPI_ON
                0242         CALL FORWARD_STEP(iloop_ocn, myTime, myIter, myThid )
                0243 #  else
                0244         myIter = nIter0 + iloop_ocn
                0245         myTime = startTime + deltaTClock *float (iloop_ocn)
85070fa4eb Jean*0246         CALL DO_THE_MODEL_IO( .FALSE., myTime, myIter, myThid )
                0247         CALL DO_WRITE_PICKUP( .FALSE., myTime, myIter, myThid )
09a6f3668a Jeff*0248 #  endif
                0249 #endif
                0250 #ifdef ML_2D
49512b0964 Jeff*0251         CALL ocean_ml(dtocn*3600.,dtatm*3600.)
09a6f3668a Jeff*0252 #endif
                0253
85070fa4eb Jean*0254       ENDDO ! ncall_ocean loop
09a6f3668a Jeff*0255
                0256 C Start of code done at the end of every coupled period
                0257
9274434acc Jean*0258 #ifdef OCEAN_3D
                0259       CALL NORM_OCN_FLUXES(myThid)
09a6f3668a Jeff*0260       CALL ATM2D_WRITE_PICKUP(.FALSE., myTime, myIter, myThid)
                0261 #endif
                0262
                0263 C
                0264 C--------------------- End of coupled period loop --------------------
                0265 C
                0266       IF (inday.EQ.0) THEN  !do this block if end-of-day
                0267 C        PRINT *,'***block at end of day'
                0268       ENDIF  !end block end-of-day
                0269
8ad2868148 Jeff*0270       IF (inyr.EQ.jdofmhr(monid)) THEN !do block if month-end
260186e531 Jeff*0271         PRINT *,'***end of month reached'
09a6f3668a Jeff*0272 #ifdef CLM
                0273 #  ifdef CPL_TEM
                0274         CALL climate2tem(monid,ndmonth(monid))
260186e531 Jeff*0275 c        PRINT *,'From driver before call tem',' idyear=',idyear
09a6f3668a Jeff*0276         CALL tem(idyear,monid-1)
                0277         CALL tem2climate(idyear,monid-1)
                0278         ch4mn=0.0
                0279         n2omn=0.0
                0280         nepmn=0.0
                0281         DO j=1,jm0
                0282           ch4mn=ch4mn+temch4(j)
                0283           n2omn=n2omn+temn2o(j)
                0284           nepmn=nepmn+temco2(j)
                0285         ENDDO
                0286 #    ifdef CPL_NEM
                0287         PRINT *,'Month=',monid
                0288         PRINT *,'CH4=',ch4mn/1.e9,' N2O=',n2omn/1.e9
49512b0964 Jeff*0289         OPEN(277,ACCESS='APPEND',FILE=fnememiss,form='unformatted'
                0290      &     ,STATUS='old')
                0291         WRITE (277) iyr,monid,ch4mn,n2omn,nepmn,
09a6f3668a Jeff*0292      &          temch4,temn2o,temco2
f6630d3a9c Jeff*0293         CLOSE(277)
09a6f3668a Jeff*0294 #    endif
                0295         DO j=1,jm0
                0296           temnep(monid,j)=temco2(j)
50a1736d54 Jean*0297         ENDDO
09a6f3668a Jeff*0298         DO j=1,jm0
                0299           antemnep(j)=antemnep(j)+temnep(monid,j)
fba33e17f0 Jeff*0300           nepan=nepan+temnep(monid,j)
09a6f3668a Jeff*0301           ch4ann=ch4ann+temch4(j)
                0302           n2oann=n2oann+temn2o(j)
50a1736d54 Jean*0303         ENDDO
09a6f3668a Jeff*0304
                0305 #  endif
                0306 #endif
                0307
                0308 #ifdef OCEAN_3D
                0309         CALL MONTH_END_DIAGS( monid, myTime, myIter, myThid)
                0310 #endif
                0311
                0312 #ifdef CPL_OCEANCO2
                0313         IF (monid.EQ.12) THEN
                0314           ocupt=ocupt*12.e-15
260186e531 Jeff*0315 C   12.e-15 from moles to Gt carbon
09a6f3668a Jeff*0316           ocuptp=ocupt
                0317           ocupt=0.0
                0318         ENDIF
                0319 #endif
                0320
                0321 #ifdef IPCC_EMI
                0322         PRINT *,'Month=',monid
                0323         nepmn=nepmn/1000.
9a05518332 Jeff*0324 C         nepmn is converted from Tg/mn to Gt/mn of C
09a6f3668a Jeff*0325         ocumn=ocumn*12.e-15
9a05518332 Jeff*0326 C         ocumn is converted from mole(C) to Gt (C)
09a6f3668a Jeff*0327 C        tnow= jyear + (jday-.5)/365.
9a05518332 Jeff*0328 C        CALL emissipcc(tnow,nepmn,ocumn,CO2,xco2ann)
260186e531 Jeff*0329          print *,nepmn,ocumn,xco2ann
9a05518332 Jeff*0330 C         ch4mn is in kg/mn of CH4
                0331 C         nepmn is in Gt/mn of C
                0332          tcumn=nepmn-ch4mn*12./16.*1.e-12
260186e531 Jeff*0333          print *,'tcumn,ocumn,xco2ann'
                0334          print *,tcumn,ocumn,xco2ann
9a05518332 Jeff*0335         CALL emissipcc_mn(tcumn,ocumn,xco2ann)
260186e531 Jeff*0336 C       CALL emissipcc_mn(nepmn,ocumn,xco2ann)
09a6f3668a Jeff*0337 #endif
                0338       ENDIF  !end block done at month-end
9274434acc Jean*0339
8ad2868148 Jeff*0340       IF (inyr .EQ. jdofmhr(12)) THEN ! do this block at year-end
09a6f3668a Jeff*0341         PRINT *,'***end of year reached'
9274434acc Jean*0342 #ifdef CPL_TEM
09a6f3668a Jeff*0343         nepan=nepan/1000.
                0344         IF (iyr.ge.1981.and.iyr.le.1990) THEN
                0345           PRINT *,'Uptake avegaging year=',iyr
                0346           nepav=nepav+nepan
                0347           aocuav=aocuav+OCUPTP
                0348           IF (iyr.eq.1990) THEN
                0349             nepav=nepav/10.
                0350             aocuav=aocuav/10.
                0351             totup=nepav+aocuav
                0352             aduptt=4.1-totup
                0353             PRINT *,' Carbon uptake for spinup'
                0354             PRINT *,' totup=',totup,' aocuav=',aocuav
                0355             PRINT *,' nepav=',nepav,' aduptt=',aduptt
                0356           ENDIF
                0357         ENDIF
                0358
                0359         IF (iyr.eq.endYear) THEN
                0360 C         NEM emissions and NEP for start of climate-chemistry run
                0361           adupt=aduptt
260186e531 Jeff*0362            CALL wr_rstrt_nem
09a6f3668a Jeff*0363         ENDIF
                0364
                0365 #endif
                0366
                0367 #ifdef ML_2D
                0368 C    Data for the restart of the 2D ML model
                0369         CALL wrrstrt_ocean
                0370 #endif
                0371
                0372 #ifdef OCEAN_3D
                0373         IF ((mod(iyr,taveDump).EQ.0).AND.(idyear.GE.taveDump)) THEN
                0374           CALL TAVE_END_DIAGS( taveDump, myTime, myIter, myThid)
                0375         ELSEIF (mod(iyr,taveDump).EQ.0) THEN
                0376           CALL TAVE_END_DIAGS( idyear, myTime, myIter, myThid)
                0377         ENDIF
8e101fde6e Jeff*0378         CALL YR_END_DIAGS(iyr,myTime,myIter,myThid)
09a6f3668a Jeff*0379 C If necessary, next line can be moved outside OCEAN_3D for IGSM2.2 cleanups
7b5b62b00f Jeff*0380         IF (iloop.EQ.nTimeSteps) CALL ATM2D_FINISH(myThid)
260186e531 Jeff*0381 #  ifdef NCEPWIND
                0382         OPEN(unit=334,file='rand_state_new.dat',status='replace')
                0383         WRITE(334,*) JSEED,IFRST,NEXTN
                0384         CLOSE(334)
                0385 #  endif
09a6f3668a Jeff*0386 #endif
                0387
d09f90e9d3 Jeff*0388 #if (defined CPL_TEM) || (defined CPL_OCEANCO2)
e72648f429 Jeff*0389         PRINT '(a6,i6,2(a5,f10.4))','Year=',iyr,
d09f90e9d3 Jeff*0390      &         ' NEP=',nepan,' OCU=',OCUPTP
                0391 #  ifdef IPCC_EMI
                0392         PRINT '(a6,i6,(a5,f10.4))','Year=',iyr,
                0393      &         ' CO2AN=',xco2ann/12.
                0394         CALL emissipcc_yr
                0395 #  endif
f345ec4dc7 Jeff*0396
d09f90e9d3 Jeff*0397 #  ifdef CPL_NEM
f345ec4dc7 Jeff*0398         PRINT *,' CH4=',ch4ann,' N2O=',n2oann
09a6f3668a Jeff*0399 #  endif
d09f90e9d3 Jeff*0400         OPEN(333,ACCESS='APPEND',FILE=caruptfile,STATUS='old')
f93afab396 Jeff*0401 #    ifndef CPL_TEM
                0402 C       For ocean carbon model only
260186e531 Jeff*0403         WRITE(333,'(i7,3e15.5)')iyr,ocuptp
                0404 #    else
                0405 #      ifndef CPL_OCEANCO2
                0406 C       For ocean TEM only
                0407         WRITE(333,'(i7,3e15.5)')iyr,nepan,nepan-1.e-12*ch4ann*12./16.
                0408 #      else
f93afab396 Jeff*0409 C        For ocean both TEM OCM
260186e531 Jeff*0410         WRITE(333,'(i7,3e15.5)')iyr,nepan,nepan-1.e-12*ch4ann*12./16.
                0411      &                  ,ocuptp
                0412 #      endif
                0413 #    endif
f6630d3a9c Jeff*0414         CLOSE(333)
d09f90e9d3 Jeff*0415 #  if (defined CPL_OCEANCO2) && (defined ML_2D)
09a6f3668a Jeff*0416         WRITE(602)iyr
                0417         CALL wrgary
                0418         CALL zerogary
d09f90e9d3 Jeff*0419 #  endif
f345ec4dc7 Jeff*0420 #endif
                0421
09a6f3668a Jeff*0422 #ifdef CPL_OCEANCO2
                0423         DO j=1,jm0
50a1736d54 Jean*0424 #  ifdef OCEAN_3D
f6630d3a9c Jeff*0425           co24ocnan(j)=co24ocnan(j)*dtatm/24.0/365.0
f345ec4dc7 Jeff*0426 #  else
49512b0964 Jeff*0427           co24ocnan(j)=co24ocnan(j)/365.0
f345ec4dc7 Jeff*0428 #  endif
09a6f3668a Jeff*0429         ENDDO
                0430         PRINT *,' CO2 for ocean model',' ncallatm=',ncall_atm
e72648f429 Jeff*0431         PRINT '(12f7.1,/,2(11f7.1,/),12f7.1)',co24ocnan
09a6f3668a Jeff*0432 #endif
f345ec4dc7 Jeff*0433
09a6f3668a Jeff*0434 #ifdef CPL_CHEM
                0435         PRINT *,' TEMUPTANN=',temuptann,' TOTAL UPTAKE='
                0436      &          ,ocuptp+temuptann
                0437 #endif
                0438       ENDIF  !year-end block
                0439
                0440       RETURN
                0441       END
260186e531 Jeff*0442
                0443 #ifdef NCEPWIND
                0444       REAL*8 FUNCTION RAND()
                0445 C
                0446 C  This function returns a pseudo-random number for each invocation.
50a1736d54 Jean*0447 C  It is a FORTRAN 77 adaptation of the "Integer Version 2" minimal
260186e531 Jeff*0448 C  standard number generator whose Pascal code appears in the article:
                0449 C
50a1736d54 Jean*0450 C     Park, Steven K. and Miller, Keith W., "Random Number Generators:
                0451 C     Good Ones are Hard to Find", Communications of the ACM,
260186e531 Jeff*0452 C     October, 1988.
                0453 C
                0454       PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773,
                0455      +           MOMDMP=2836)
                0456 C
                0457       COMMON  /SEED/JSEED,IFRST,NEXTN
                0458       INTEGER JSEED,IFRST,NEXTN
                0459       INTEGER HVLUE, LVLUE, TESTV
                0460 C
                0461       IF (IFRST .EQ. 0) THEN
                0462         NEXTN = JSEED
                0463         IFRST = 1
                0464       ENDIF
                0465 C
                0466       HVLUE = NEXTN / MOBYMP
                0467       LVLUE = MOD(NEXTN, MOBYMP)
                0468       TESTV = MPLIER*LVLUE - MOMDMP*HVLUE
                0469       IF (TESTV .GT. 0) THEN
                0470         NEXTN = TESTV
                0471       ELSE
                0472         NEXTN = TESTV + MODLUS
                0473       ENDIF
                0474       RAND = REAL(NEXTN)/REAL(MODLUS)
                0475 C
                0476       RETURN
                0477       END
                0478 #endif