mo_nudging.f90
MODULE mo_nudging
!======================================================================
!
5: ! J. Feichter Uni Hamburg Jan 91 and Apr 93
! W. May DMI-Copenhagen Mar 97
! M. Stendel MPI-Hamburg Sep 96 and May 97
! H.-S. Bauer MPI-Hamburg Jul 98
! I. Kirchner MPI-Hamburg Nov 98 and Oct 99
10: !
! ****************** Interface to ECHAM4 ******************************
!
! *ECHAM* *NUDGING*
!
15: ! INICTL ---> NudgingInit(0) initialize nudging
! CONTROL --> NudgingInit(10)
! NudgingInit(1)
!
! STEPON ---> NudgingReadSST read new global sst field
20: ! |
! +-> Nudging perform nudging
! | +--->GetNudgeData read nudging data sets
! | +-->ReadOneBlock read one data time step
! | +-->OpenOneBlock open new data block
25: ! | +-->CloseBlock close data block
! |
! +-> NudgingOut store nudging terms (accu+inst)
! | +-->NudgingStoreSP store one spectral array
! |
30: ! +-> NudgingRerun(-1) read/write nudging buffer
!
! GPC ------> NudgingSSTnew reset sst at latitude circle
!
!======================================================================
35: !
! This version is for real time. SST should be read once a day
! (to be set in READSST).
! --------------------------------------------------------------------
! Following Krishnamurti et al. (1991), Tellus 43AB, 53-81.
40: !
! Now nudge the coefficients and compute the tendencies.
! No nudging of first coefficient, since we do not want to force the
! global mean toward an unbalanced state.
!
45: ! This is "pure" Krishnamurti, which means that there are no other
! relaxation terms on rhs!
!
! (Eq. 7.3): A(t+dt) = (A (t+dt) +2dt*N * A (t+dt)) / (1+2dt*N),
! * 0
50: !
! where A (t+dt) is a predicted value of A (t+dt) prior to nudging.
! *
! A (t+dt) is a future value which the nudging is aimed at.
! 0 (in the twin case that is the observed value
55: ! at the later time)
! N is nudging coefficient
!
! the method is changed
!
60: ! NEW = A*OLD + B*OBS
!
! implicit --> A = 1/(1-2*dt*N) B = 2*dt*N/(1-2*dt*N)
!
! explicit --> A = 1 - 2*dt*N B = 2*dt*N
65: !
! --------------------------------------------------------------------
USE mo_kind, ONLY: cp
USE mo_exception, ONLY: finish
70: USE mo_doctor, ONLY: nout, nin
USE mo_start_dataset, ONLY: lres, nstart, nstep
USE mo_nudging_buffer
USE mo_memory_g3a, ONLY: auxil1m, auxil2m, tsm, tsm1m, slmm
USE mo_memory_sp, ONLY: sd, svo, stp
75: USE mo_nmi, ONLY: NMI_Make
USE mo_year, ONLY: iymd2c, cd2dat, ic2ymd, isec2hms, im2day
USE mo_physc2, ONLY: ctfreez
USE mo_constants, ONLY: dayl
USE mo_truncation, ONLY: nmp, nnp
80: USE mo_control, ONLY: nlon, nlev, nsp, nlevp1, ncbase, nm, &
& ntbase, dtime, twodt, nresum, ngl, lnudge, ltdiag, nmp1, nn, &
& ncdata, ntdata, lnmi, nrow, n2sp, nptime, vct
USE mo_post, ONLY: nunitsp, nunitdf, lppt, lppp, lppvo, lppd
85: USE mo_grib, ONLY: kgrib, kleng, klengo, ksec0, ksec1, ksec2_sp, psec2, ksec3, psec3, &
& ksec4, kword, nbit, iobyte, nudging_table, year, month, day, hour, minute, century, &
& code_parameter, level_type, level_p2, set_output_time
IMPLICIT NONE
90:
! external nudging weigths
INTEGER, PARAMETER :: nmaxc = 50
REAL :: nudgd(nmaxc), nudgv(nmaxc), nudgt(nmaxc), nudgp
REAL, PARAMETER :: nfact=1.e-5 ! scaling factor for nudging coefficients
95: ! internal weights for observed data
REAL, POINTER, DIMENSION(:) :: nudgdo, nudgvo, nudgto
REAL :: nudgpo
! local fields with nudging coefficients
REAL, POINTER, DIMENSION(:,:,:) :: nudgda, nudgdb, nudgva, nudgvb, nudgta, nudgtb
100:
! true = use first date from nudging data block
! false = use information from history data (default)
LOGICAL :: lnudgini
105: ! true = implicit nudging (default)
! false = explicit nudging
LOGICAL :: lnudgimp
! nudgsmin -1 .......... incl. global mean
110: ! 0 ... nmp1-1 zonal wavenumber
INTEGER :: nudgsmin, nudgsmax, nudgmin, nudgmax, nudglmin, nudglmax
REAL :: nudgdamp
! wave specific truncation
115: ! 0 ... use all coefficients of active wave no.
! 1 ... cut off meridional index higher than wave no. limit all waves
! 2 ... cut off except wave 0
INTEGER :: nudgtrun
INTEGER, POINTER, DIMENSION(:,:,:) :: &
120: & flagn, & ! set to 1 inside the nudging region
& flago ! set to 1 outside nudging region
! twin data processing
LOGICAL :: ltwin
125: INTEGER :: ntwinc, ntwoff
! sst data usage
INTEGER :: nsstinc, nsstoff
130: ! control the output
LOGICAL :: lnudgwobs
INTEGER :: ndunit(13) ! div, vor, temP time1 == 1 2 3
! time2 == 4 5 6
! time3 == 7 8 9
135: ! sst t1 t2 t3 == 10 11 12
! for restart file == 13
! internal file pointer of nudging data blocks and sst fields
! DIV, VOR, TEM+P, SST, RERUN
140: INTEGER (cp) :: kfiled, kfilev, kfilet, kfiles, kfiler
CHARACTER*10 :: cfile
INTEGER :: ndgblock ! block pointer nudging data
INTEGER :: sstblock ! number of sst file
145: ! time/date marker
INTEGER :: ndgstep1
INTEGER :: ihead1d, ihead1t ! YYYYMMDD HHmmss
INTEGER :: ndgstep2
INTEGER :: ihead2d, ihead2t ! YYYYMMDD HHmmss
150:
LOGICAL :: lsstn ! true if new sst neccessary
INTEGER :: iheads ! YYYYMMDDHH sst time step
REAL, POINTER, DIMENSION (:,:) :: sstn ! (nlon,ngl)
155: ! code numbers : Pins,Tins,Dins,Vins,Pacc,Tacc,Dacc,Vacc
! this is for GRIB1 special code page usage
! INTEGER, PARAMETER :: ndgcode(20)=(/129, 130, 131, 132, 133, 134, 135, 136, &
!& 137, 138, 139, 140, 141, 142, 143, 144,&
!& 145, 146, 147, 148/)
160: INTEGER, PARAMETER :: ndgcode(20)=(/111, 112, 113, 114, 115, 116, 117, 118, &
& 21, 22, 23, 24, 25, 26, 27, 28, &
& 31, 32, 33, 34/)
! first set: nudging term
! second set: difference to observations outside the nudging range
165: ! third set: observations instantaneous
! nudging tendencies (instantaneous and accumulated)
REAL, POINTER, DIMENSION(:,:,:) :: sdten, svten, stten, sdtac, svtac, sttac
170: ! write more messages
LOGICAL :: lnudgdbx
CONTAINS
!======================================================================
175:
SUBROUTINE NudgingInit(itype)
! setup nudging parameter and fields
USE mo_mpi
180:
INCLUDE 'ndgctl.inc'
INTEGER :: i, kret, iday, id, im, iy, isec, ihms, iymd, jk
INTEGER :: ihead(8), itype, lday, is1, is2, jm, ii, sum1, sum2, sum3
185: CHARACTER (8) :: yhead(8)
! Intrinsic functions
INTRINSIC MOD, MAX, MIN
190: WRITE(nout,*) 'NudgingInit:'
IF (lnudge) THEN
!-- A. nudging
195: SELECT CASE(itype)
!-- 1. namelist setup and time adjustment
CASE(0) ! read namelist and set time syncronization
200:
WRITE(nout,*) 'Nudging Version Rev. 2.3 10-OCT-99 (kirchner@dkrz.de)'
!----------------------------------------------------------------
! preset variables
205: lnudgdbx = .FALSE.
nsstinc = 24 ! sst data given at every 24 hour
nsstoff = 12 ! 12 UTC sst data used
nudgsmin = 0 ! skip global average
nudgsmax = nn ! all other spectral components used
210: nudglmin = 1 ! from top
nudglmax = 19 ! to bottom
nudgdamp = 1.0 ! no weigthing between two time steps
nudgtrun = 0 ! use all zonal indexes
215: ! set file units
! first nudging data block
! second nudging data block
! third nudging data block
DO i=1,9
220: ndunit(i) = 80+i
END DO
! sst data block
ndunit(10) = 46
ndunit(11) = 47
225: ndunit(12) = 48
! restart file
ndunit(13) = 40
! set nudging coefficients to standard (Jeuken et al. 1996)
230: IF (nlev > nmaxc) CALL finish('NudgingInit','number of levels too large')
nudgd(:) = 5.0
nudgv(:) = 10.0
nudgt(:) = 1.0
nudgp = 10.0
235:
! allocate weights of observations
ALLOCATE(nudgdo(nlev)); nudgdo(:) = 0.0
ALLOCATE(nudgvo(nlev)); nudgvo(:) = 0.0
ALLOCATE(nudgto(nlev)); nudgto(:) = 0.0
240: nudgpo = 0.0
lnudgini = .FALSE. ! use rerun data information
lnudgimp = .TRUE. ! use implicit nudging methode
lnudgwobs = .FALSE. ! no additional output
245: ltwin = .FALSE. ! read nudging data as default
ntwinc = 3 ! twin-data output in 3 hour intervalls
ntwoff = 0 ! twin output start at 00 UTC
!----------------------------------------------------------------
250: ! read nudging namelist
READ (nin,ndgctl)
!----------------------------------------------------------------
255: ! correct nudging namelist parameters
IF (.NOT.lres) lnudgini = .TRUE.
IF (ltwin) lnudgini = .FALSE.
260: ntwinc = MAX(MIN(ntwinc,23),1)
ntwoff = MAX(MIN(ntwoff,23),0)
WRITE(nout,*) ' switches: ', &
& ' LNUDGINI = ',lnudgini, &
265: & ' LNUDGIMP = ',lnudgimp, &
& ' LNUDGDBX = ',lnudgdbx, &
& ' LNUDGWOBS = ',lnudgwobs
nudgdamp = MIN(MAX(nudgdamp,0.0),1.0)
270:
WRITE(nout,*) ' maximal damping in the middle of data time steps ',nudgdamp
WRITE(nout,*) ''
WRITE(nout,*) ' NUDG data block 1 from units ',(ndunit(i),i= 1, 3)
WRITE(nout,*) ' NUDG data block 2 from units ',(ndunit(i),i= 4, 6)
275: WRITE(nout,*) ' NUDG data block 3 from units ',(ndunit(i),i= 7, 9)
WRITE(nout,*) ' NUDG read sst from units ',(ndunit(i),i=10,12)
WRITE(nout,*) ' NUDG restart file unit ', ndunit(13)
IF (ltwin) THEN
280: WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(11)
CALL pbopen(kfiles,cfile,'w',kret)
WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(4)
CALL pbopen(kfiled,cfile,'w',kret)
285:
WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(5)
CALL pbopen(kfilev,cfile,'w',kret)
WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(6)
290: CALL pbopen(kfilet,cfile,'w',kret)
WRITE(nout,*) ' TWIN mode use block 2 for data storage'
WRITE(nout,*) ' first time at day ',ntwoff,' UTC increment ',ntwinc,' hours'
nsstinc = ntwinc
295: nsstoff = ntwoff
ELSE
ntwinc = -1
ntwoff = -1
END IF
300:
!----------------------------------------------------------------
! synchronize date and time
sstblock = -1 ! start in first sst file
305:
IF (lnudgini) THEN
! start with old rerun data or initial data
! adjust date to first record block 2, it should be 00UTC
310:
WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(4)
ndgblock = 2
! read first time step of nudging data block 2
315: !
CALL pbopen(kfiled,cfile,'r',kret)
CALL pbread(kfiled,yhead(1),8*8,kret)
CALL util_i8toi4 (yhead(1), ihead(1),8)
CALL pbclose(kfiled,kret)
320:
! set one time step before 00UTC
! the initial time offset is (nstep+1)
IF (lres) THEN
! use old rerun data
325: nstep = 1
ELSE
! use initial data
nstep = 0
ENDIF
330:
! calculate new initial date and time, ncbase and ntbase
!
! YYYYMMDD ihead(3), HH ihead(4)
!
335: ncbase = iymd2c(ihead(3)) ! set initial day in the century
ntbase = ihead(4)*3600 - (nstep+1)*dtime
IF (ntbase < 0) THEN
! correct time at midnight
ncbase = ncbase - 1
340: ntbase = dayl + ntbase
ENDIF
! reset verification time
ncdata = ncbase
ntdata = ntbase
345: !
CALL cd2dat(ncbase,id,im,iy)
WRITE(nout,*) ' reset NCBASE= ',ncbase,' NTBASE= ',ntbase,' NSTEP= ',nstep
WRITE(nout,*) ' corresponding date Year ',iy,' Month ',im,' Day',id
!
350: ! set first date of nudging data
iday = ncbase + (ntbase+dtime*(nstep+1))/dayl + 0.01
CALL cd2dat(iday,id,im,iy)
WRITE(nout,*) ' first nudging date Year ',iy,' Month ',im,' Day',id,' 00UTC'
355:
nresum = nstep
nstart = 0
ELSE IF (.NOT.ltwin) THEN
360:
! read all nudging data blocks
ndgblock = 1
! set first date of nudging data
365: iday = ncbase + (ntbase+dtime*(nstep+1))/dayl + 0.01
iymd = ic2ymd(iday)
isec = INT(MOD((ntbase + dtime*(nstep+1)),dayl))
ihms = isec2hms(isec)
370: !
WRITE(nout,*) ' rerun NCBASE= ',ncbase,' NTBASE= ',ntbase,' NSTEP= ',nstep
WRITE(nout,*) ' first nudging date expected ... YYMMDD ' &
& ,iymd,' HHMMSS = ',ihms
375: WRITE(nout,*) ' use time information from restart files'
ENDIF
! reset nudging data block pointer and time information
380: ihead1d = 0
ihead1t = 0
ihead2d = 0
ihead2t = 0
385: ! initialize sst field pointer
nsstinc = MAX(nsstinc,0)
IF (nsstinc==0) THEN
WRITE(nout,*) ' WARNING :: NSSTINC is zero, use standard ECHAM SST'
iheads = -99
390: ELSE
ALLOCATE (sstn(nlon,ngl))
iheads = 0
ENDIF
395: WRITE(nout,*) ' Nudging ... basic initialization done'
!-- 2. Setup nudging coefficients and work space
CASE(1)
400:
! correct spectral limits
nudgsmin = MIN(MAX(nudgsmin,-1),nm)
nudgsmax = MIN(MAX(nudgsmax, 0),nm)
nudgsmax = MAX(nudgsmax,nudgsmin)
405:
nudgtrun = MAX(MIN(nudgtrun,2),0)
IF (ltwin) THEN
nudgsmin = -1
410: nudgsmax = nm
nudgtrun = 3
WRITE (nout,*) ' twin mode on'
ELSE
415:
WRITE(nout,*) ' nudge zonal wavenumers (including): ',nudgsmin,' ... ',nudgsmax
WRITE(nout,*) ' use truncation type ',nudgtrun
! set spectral index of nudging region
420: IF (nudgsmin<0) THEN
WRITE(nout,*) ' nudge also global average'
nudgmin = 1
nudgsmin = 0
ELSE IF (nudgsmin == 0) THEN
425: ! nudge wave number 0 without global average
nudgmin = 2
ELSE
nudgmin = nmp(nudgsmin+1)+1
ENDIF
430: nudgmax = nmp(nudgsmax+1)+nnp(nudgsmax+1)
nudgmax = MAX(nudgmax,nudgmin)
! set vertical region
nudglmin = MIN(MAX(nudglmin,1),nlev)
435: nudglmax = MIN(MAX(nudglmax,1),nlev)
nudglmax = MAX(nudglmax,nudglmin)
! set region flags
ALLOCATE (flagn(nlevp1,2,nsp))
440: ALLOCATE (flago(nlevp1,2,nsp))
flagn(:,:,:) = 0
flago(:,:,:) = 1
SELECT CASE(nudgtrun)
445: CASE (0)
WRITE(nout,*) ' truncation ... use all zonal numbers'
flagn(nudglmin:nudglmax,:,nudgmin:nudgmax) = 1
flagn(nlevp1 ,:,nudgmin:nudgmax) = 1
flago(nudglmin:nudglmax,:,nudgmin:nudgmax) = 0
450: flago(nlevp1 ,:,nudgmin:nudgmax) = 0
CASE (1)
WRITE(nout,*) ' truncation ... cut off zonal numbers'
DO jm=nudgsmin,nudgsmax
is1 = nmp(jm+1) + 1
455: is2 = nmp(jm+1) + 1 + nudgsmax - jm
flagn(nudglmin:nudglmax,:,is1:is2) = 1
flagn(nlevp1 ,:,is1:is2) = 1
flago(nudglmin:nudglmax,:,is1:is2) = 0
flago(nlevp1 ,:,is1:is2) = 0
460: END DO
CASE (2)
WRITE(nout,*) ' truncation ... cut off except wave 0'
DO jm=nudgsmin,nudgsmax
IF (jm==0) THEN
465: is1 = nmp(1) + 1
is2 = nmp(1) + nnp(1)
ELSE
is1 = nmp(jm+1) + 1
is2 = nmp(jm+1) + 1 + nudgsmax - jm
470: END IF
flagn(nudglmin:nudglmax,:,is1:is2) = 1
flagn(nlevp1 ,:,is1:is2) = 1
flago(nudglmin:nudglmax,:,is1:is2) = 0
flago(nlevp1 ,:,is1:is2) = 0
475: END DO
END SELECT
! check flag field
sum1 = 0
480: sum2 = 0
sum3 = 0
DO i=1,nlevp1
DO ii=1,nsp
sum1 = sum1 + flagn(i,1,ii)+flagn(i,2,ii)
485: sum2 = sum2 + flago(i,1,ii)+flago(i,2,ii)
sum3 = sum3 + flagn(i,1,ii)*flago(i,1,ii)+flagn(i,2,ii)*flago(i,2,ii)
END DO
END DO
IF ((sum1+sum2)/=(n2sp*nlevp1)) THEN
490: WRITE(nout,*) ' WARNING nudging flag field not correct'
WRITE(nout,*) ' FLAGN = ',sum1,' FLAGO = ',sum2,' O*N = ',sum3,' NO = ',n2sp*nlevp1
END IF
! rescale nudging coefficients and set level window
495: IF (nudglmin > 1) THEN
nudgd(1:nudglmin-1) = 0.0
nudgv(1:nudglmin-1) = 0.0
nudgt(1:nudglmin-1) = 0.0
ENDIF
500: IF (nudglmax < nlev) THEN
nudgd(nudglmax+1:nlev) = 0.0
nudgv(nudglmax+1:nlev) = 0.0
nudgt(nudglmax+1:nlev) = 0.0
ENDIF
505:
!----------------------------------------------------------------
! prepare nudging weights level dependend
!
! set weights for explicit nudging
510: ! new = N(read)*2*dt*nfact
! rescale coefficients
nudgdo(:) = nudgd(1:nlev) * nfact * twodt
nudgvo(:) = nudgv(1:nlev) * nfact * twodt
nudgto(:) = nudgt(1:nlev) * nfact * twodt
515: nudgpo = nudgp * nfact * twodt
! set weights for model data
nudgd(1:nlev) = 1.-nudgdo(:)
nudgv(1:nlev) = 1.-nudgvo(:)
520: nudgt(1:nlev) = 1.-nudgto(:)
nudgp = 1.-nudgpo
! convert coefficients for implicit nudging
IF (lnudgimp) THEN
525: ! using implicit nudging reset weights
nudgd(1:nlev) = 1./(1.+nudgdo(:))
nudgv(1:nlev) = 1./(1.+nudgvo(:))
nudgt(1:nlev) = 1./(1.+nudgto(:))
nudgp = 1./(1.+nudgpo)
530:
! reset weights for observations
nudgdo(:) = nudgdo(:) * nudgd(1:nlev)
nudgvo(:) = nudgvo(:) * nudgv(1:nlev)
nudgto(:) = nudgto(:) * nudgt(1:nlev)
535: nudgpo = nudgpo * nudgp
ENDIF
WRITE(nout,*) &
& ' nudging coefficients (1/sec) and weights for model and observations (%)'
540: WRITE(nout,'(a5,3a25,7x)') 'lev','DIV','VOR','TEM'
IF (lnudgimp) THEN
WRITE(nout,*) ' use IMPLICIT nudging'
545: DO i=1,nlev
WRITE(nout,'(i5,3(e12.3,2f10.3))') i,&
& (1.0-nudgd(i))/(nudgd(i)*twodt),nudgd(i)*100.0,nudgdo(i)*100.0,&
& (1.0-nudgv(i))/(nudgv(i)*twodt),nudgv(i)*100.0,nudgvo(i)*100.0,&
& (1.0-nudgt(i))/(nudgt(i)*twodt),nudgt(i)*100.0,nudgto(i)*100.0
550: ENDDO
WRITE(nout,'(a,e12.3,2f10.3)') ' pressure = ',&
& (1.0-nudgp)/(nudgp*twodt),nudgp*100.0,nudgpo*100.0
ELSE
555:
WRITE(nout,*) ' use EXPLICIT nudging'
DO i=1,nlev
WRITE(nout,'(i5,3(e12.3,2f10.3))') i,&
& nudgdo(i)/twodt,nudgd(i)*100.0,nudgdo(i)*100.0,&
560: & nudgvo(i)/twodt,nudgv(i)*100.0,nudgvo(i)*100.0,&
& nudgto(i)/twodt,nudgt(i)*100.0,nudgto(i)*100.0
ENDDO
WRITE(nout,'(a,e12.3,2f10.3)') ' pressure = ',&
& nudgpo/twodt,nudgp*100.0,nudgpo*100.0
565:
ENDIF
WRITE(nout,*) ' nudging between (including) LEVEL = ',nudglmin,' and ',nudglmax
WRITE(nout,*) ' nudging between (including) SPECC = ',nudgmin, ' and ',nudgmax
570:
!----------------------------------------------------------------
! set work space nudging coefficients
!
ALLOCATE(nudgda(nlev ,2,nsp)); nudgda(:,:,:) = 0.0
575: ALLOCATE(nudgdb(nlev ,2,nsp))
ALLOCATE(nudgva(nlev ,2,nsp)); nudgva(:,:,:) = 0.0
ALLOCATE(nudgvb(nlev ,2,nsp))
ALLOCATE(nudgta(nlevp1,2,nsp)); nudgta(:,:,:) = 0.0
ALLOCATE(nudgtb(nlevp1,2,nsp))
580:
! set first part of coefficients
! second part set in Nudging
DO jk=1,nlev
nudgdb(jk,:,:) = nudgdo(jk)*flagn(jk,:,:)
585: nudgvb(jk,:,:) = nudgvo(jk)*flagn(jk,:,:)
nudgtb(jk,:,:) = nudgto(jk)*flagn(jk,:,:)
ENDDO
nudgtb(nlevp1,:,:) = nudgpo*flagn(nlevp1,:,:)
590: WRITE(nout,*) ''
!----------------------------------------------------------------
!
IF (lnudgini) THEN
595: ! time step correction
! the initial time offset is (nstep+1)
IF (lres) THEN ! use old rerun data
nstep = 1
ELSE ! use initial data
600: nstep = 0
ENDIF
nresum = nstep
nstart = 0
ELSE
605: ! reset the input data block pointer
iday = ncbase + (ntbase+dtime*(nstep+1))/dayl + 0.01
CALL cd2dat(iday,id,im,iy)
lday = im2day(im,iy)
IF (id<lday) THEN
610: ndgblock = 2
WRITE(nout,*) 'NudgingInit: skip first nudging data block'
ENDIF
ENDIF
615: !----------------------------------------------------------------
! initialize nudging fields
ALLOCATE (sdobs1(nlev ,2,nsp)); sdobs1(:,:,:) = 0.
ALLOCATE (svobs1(nlev ,2,nsp)); sdobs1(:,:,:) = 0.
620: ALLOCATE (stobs1(nlevp1,2,nsp)); sdobs1(:,:,:) = 0.
ALLOCATE (sdobs2(nlev ,2,nsp)); sdobs2(:,:,:) = 0.
ALLOCATE (svobs2(nlev ,2,nsp)); sdobs2(:,:,:) = 0.
ALLOCATE (stobs2(nlevp1,2,nsp)); sdobs2(:,:,:) = 0.
625:
ALLOCATE (sdobs(nlev ,2,nsp)); sdobs(:,:,:) = 0.
ALLOCATE (svobs(nlev ,2,nsp)); sdobs(:,:,:) = 0.
ALLOCATE (stobs(nlevp1,2,nsp)); sdobs(:,:,:) = 0.
630: ALLOCATE (sdten(nlev ,2,nsp))
ALLOCATE (svten(nlev ,2,nsp))
ALLOCATE (stten(nlevp1,2,nsp))
ALLOCATE (sdtac(nlev ,2,nsp)); sdtac(:,:,:) = 0.0
635: ALLOCATE (svtac(nlev ,2,nsp)); svtac(:,:,:) = 0.0
ALLOCATE (sttac(nlevp1,2,nsp)); sttac(:,:,:) = 0.0
WRITE (nout,*) 'NudgingInit: memory allocated'
640: ! reload old restart data
IF(.NOT.lnudgini) CALL NudgingRerun(1)
END IF
645: !-- 3. reset NSTEP after RESTART
CASE(10)
IF (lnudgini.AND.(.NOT.lres)) THEN
650: nstep = 1
nresum = nstep
ENDIF
END SELECT
655:
ELSE
!-- B. no nudging
660: !----------------------------------------------------------------
! read nudging namelist
IF (p_parallel) THEN
IF (p_parallel_io) THEN
WRITE(nout,*) ' skip namelist'
665: READ (nin,ndgctl)
ENDIF
ELSE
WRITE(nout,*) ' skip namelist'
READ (nin,ndgctl)
670: ENDIF
ENDIF
RETURN
675:
END SUBROUTINE NudgingInit
!======================================================================
SUBROUTINE Nudging
680:
! -------------------------------------------------------------------
! New 'nudging' routine for ECHAM4 to adjust the following
! meteorological variables to observations by means of
! Newtonian relaxation: - divergence
685: ! - vorticity
! - temperature
! - log surface pressure
!
! -------------------------------------------------------------------
690: ! *NUDGING* IS CALLED FROM *STEPON*
! -------------------------------------------------------------------
REAL :: weighto, weightm, tfact1, rtwodt, tfact2
INTEGER :: jk
695:
#ifdef sun
INTEGER :: it0, it1
! External functions
EXTERNAL clockm
700: #endif
#ifdef sun
IF (lnudgdbx) CALL clockm(it0)
#endif
705:
IF (ltwin) THEN
CALL WriteTwin
710: ELSE
! initalize instantaneous tendencies
sdten(:,:,:) = 0.0
svten(:,:,:) = 0.0
715: stten(:,:,:) = 0.0
rtwodt = 1./twodt
! Read new data block
720: CALL GetNudgeData
! Interpolate the data
!
! --------------------------------------------------------------------
725: ! Read new date from analyses if necessary. For twin experiments
! in T42 (i.e. timestep=24 min) this is the case every 30 timesteps
! (24min*30=720min=12 hours).
!
! now independend of available time stepping
730: !
! --------------------------------------------------------------------
! linear time interpolation
tfact1 = REAL(nstep + 1 - ndgstep1)/REAL(ndgstep2 - ndgstep1)
735: tfact2 = 1.-tfact1
! time damping weight
weighto = nudgdamp + (1.-nudgdamp) * &
& ABS( REAL(ndgstep1+ndgstep2-2*(nstep+1)) )/REAL(ndgstep2-ndgstep1)
740: weightm = 1.-weighto
! --------------------------------------------------------------------
! Linear interpolation between the intervals of analyses
! --------------------------------------------------------------------
745:
sdobs(:,:,:) = sdobs1(:,:,:) * tfact2 + sdobs2(:,:,:) * tfact1
svobs(:,:,:) = svobs1(:,:,:) * tfact2 + svobs2(:,:,:) * tfact1
stobs(:,:,:) = stobs1(:,:,:) * tfact2 + stobs2(:,:,:) * tfact1
750: IF (lnudgdbx) WRITE(nout,*) &
& 'Nudging: time weights O= ',weighto,' M= ',weightm, &
& ' time factors T1= ',tfact1,' T2= ',tfact2,' at NSTEP+1 = ',nstep+1
IF (lnmi) CALL NMI_Make(3) ! prepare NMI filtered field
755:
! set coefficients
DO jk=1,nlev
nudgda(jk,:,:) = nudgd(jk)
nudgva(jk,:,:) = nudgv(jk)
760: nudgta(jk,:,:) = nudgt(jk)
END DO
nudgta(nlevp1,:,:) = nudgp
nudgda(:,:,:) = flago(1:nlev,:,:) + nudgda(:,:,:)*flagn(1:nlev,:,:) + weightm*nudgdb(:,:,:)
765: nudgva(:,:,:) = flago(1:nlev,:,:) + nudgva(:,:,:)*flagn(1:nlev,:,:) + weightm*nudgvb(:,:,:)
nudgta(:,:,:) = flago( : ,:,:) + nudgta(:,:,:)*flagn( : ,:,:) + weightm*nudgtb(:,:,:)
! --------------------------------------------------------------------
! Set the complex part of the global means to 0.
770: ! --------------------------------------------------------------------
sdobs(:,2,1) = 0.
svobs(:,2,1) = 0.
stobs(:,2,1) = 0.
775: !--------------------------------------------------------------------
! Now perform Newtonian relaxation.
!
! Apply full nudging for all scales. For T106 run, it may be
! necessary to decrease weight for smaller scales.
780: ! Originally, there was full nudging up to T42 (NSP<=946), 50%
! at T63 (NSP=2080) and no nudging for NSP>2628 (T71).
!--------------------------------------------------------------------
!
! --------------------------------------------------------------------
785: ! Tendencies are defined as the difference between the
! nudged and the original fields
!
! T_tendency := T_nudging - T_original, ergo
! >>> T_nudging = T_original + T_tendency <<<
790: ! --------------------------------------------------------------------
! store old model values
sdten(:,:,:) = sd (:,:,:)
svten(:,:,:) = svo(:,:,:)
795: stten(:,:,:) = stp(:,:,:)
! calculate new model value as linear combination of old model value and observations
! calculate in nudging window (nudglmin:nudglmax,nudgmin:nudgmax)
sd(:,:,:) = nudgda(:,:,:)*sd (:,:,:) + weighto*nudgdb(:,:,:)*sdobs(:,:,:)
800: svo(:,:,:) = nudgva(:,:,:)*svo(:,:,:) + weighto*nudgvb(:,:,:)*svobs(:,:,:)
stp(:,:,:) = nudgta(:,:,:)*stp(:,:,:) + weighto*nudgtb(:,:,:)*stobs(:,:,:)
! compose diagnostics with flag-fields
!
805: ! INSIDE
! store nudging term , tendency term, (NEW-OLD)/2DT in UNIT/sec
!
! OUTSIDE
! calculate difference outside the nudging range, absolut value
810: ! MODEL - OBSERVED in UNIT
!
sdten(:,:,:) = &
& flagn(1:nlev,:,:)*(sd (:,:,:) - sdten(:,:,:))*rtwodt + &
& flago(1:nlev,:,:)*(sd (:,:,:) - sdobs(:,:,:))
815:
svten(:,:,:) = &
& flagn(1:nlev,:,:)*(svo(:,:,:) - svten(:,:,:))*rtwodt + &
& flago(1:nlev,:,:)*(svo(:,:,:) - svobs(:,:,:))
820: stten(:,:,:) = &
& flagn(:,:,:)*(stp(:,:,:) - stten(:,:,:))*rtwodt + &
& flago(:,:,:)*(stp(:,:,:) - stobs(:,:,:))
! accumulation of diagnostics
825:
sdtac(:,:,:) = sdtac(:,:,:) + sdten(:,:,:)
svtac(:,:,:) = svtac(:,:,:) + svten(:,:,:)
sttac(:,:,:) = sttac(:,:,:) + stten(:,:,:)
830: ENDIF
#ifdef sun
IF (lnudgdbx) THEN
CALL clockm(it1)
835: WRITE (nout,*) 'Nudging: performance at ', nstep, it1-it0,' msec'
END IF
#endif
RETURN
840:
END SUBROUTINE Nudging
!======================================================================
SUBROUTINE GetNudgeData
845:
! read new data for nudging, dependend on the time step
INTEGER :: iday, sheadd, sheadt, isec, itime
850: ! Intrinsic functions
INTRINSIC MOD
! set expected header
itime = ntbase + dtime*(nstep+1)
855:
iday = ncbase + itime/dayl + 0.01
sheadd = ic2ymd(iday) ! get year, months and day
isec = MOD(itime,INT(dayl))
860: sheadt = isec2hms(isec) ! get hour, minute and second
IF (lnudgdbx) WRITE(nout,*) &
& 'GetNudgeData: expect date ',sheadd,' time ',sheadt
865: IF (ihead1d == 0) THEN
! after initialization load first record
CALL OpenOneBlock
CALL ReadOneBlock
IF (lnudgini) THEN
870: IF ((ihead2d/=sheadd).OR.(ihead2t/=sheadt)) THEN
!
! first record in block 2 don't correspond to shead
CALL finish('GetNudgeData','nuding data error in block 2')
ENDIF
875: ELSE
! search next possible nudging data set larger than the given date
!
! first record is larger than expected
IF ((ihead2d>sheadd).OR.((ihead2d==sheadd).AND.(ihead2t>sheadt))) &
880: & CALL finish('GetNudgeData','nudging data DATE TIME error')
! found date larger than expected date
DO
! read until date is larger than expected date
885: ihead1d = ihead2d
ihead1t = ihead2t
ndgstep1 = ndgstep2
sdobs1(:,:,:) = sdobs2(:,:,:)
svobs1(:,:,:) = svobs2(:,:,:)
890: stobs1(:,:,:) = stobs2(:,:,:)
CALL ReadOneBlock
IF ((ihead2d>sheadd).OR.((ihead2d==sheadd).AND.(ihead2t>sheadt))) EXIT
ENDDO
895: ENDIF
ENDIF
IF ((ihead2d<sheadd).OR.((ihead2d==sheadd).AND.(ihead2t<=sheadt))) THEN
!
900: ! actual time step is outside the availabe nudging data range
! read next nudging data set
ihead1d = ihead2d
ihead1t = ihead2t
ndgstep1 = ndgstep2
905: sdobs1(:,:,:) = sdobs2(:,:,:)
svobs1(:,:,:) = svobs2(:,:,:)
stobs1(:,:,:) = stobs2(:,:,:)
CALL ReadOneBlock
910: IF (lnudgdbx) WRITE(nout,*) &
& ' interpolate between ',ihead1d,ihead1t,ndgstep1 &
& ,' and ',ihead2d,ihead2t,ndgstep2
ENDIF
915:
RETURN
END SUBROUTINE GetNudgeData
!======================================================================
920:
SUBROUTINE ReadOneBlock
! read next nudging data time step
925: INTEGER :: kret, ilens1, ilens2, iilen, iday
INTEGER :: iheadd(8), iheadv(8), iheadt(8)
CHARACTER (8) :: yhead(8)
REAL, ALLOCATABLE :: phbuf(:), zhbuf(:)
930: ! Intrinsic functions
INTRINSIC MAX, RESHAPE
#ifdef CRAY
INTEGER :: iamount
935: EXTERNAL util_reshape
#endif
! --------------------------------------------------------------------
! Read analyses or model data (for twin experiment), respectively.
940: ! Here the READ statement has to be changed compared to the original
! version, since data fields are written with just one WRITE
! statement in "PREPARE_ECHAM" instead of one for each level.
! --------------------------------------------------------------------
945: ilens1 = n2sp*nlev
ilens2 = n2sp*nlevp1
iilen = MAX(ilens1,ilens2)
ALLOCATE(phbuf(iilen))
950: #ifndef CRAY
ALLOCATE(zhbuf(iilen))
zhbuf(:) = 0.0
#endif
955: DO
! read next nudging data set header
CALL pbread(kfiled,yhead(1),8*8,kret)
IF (kret/=8*8) THEN
! end of actual block reached
960: CALL CloseBlock
! open next data block
ndgblock = ndgblock + 1
CALL OpenOneBlock
CYCLE ! goto the beginning of the do-loop
965: ENDIF
CALL util_i8toi4 (yhead(1), iheadd(1),8)
CALL pbread(kfilev,yhead(1),8*8,kret)
IF (kret/=8*8) CALL finish('ReadOneBlock','nudging data inconsistent VOR')
970: CALL util_i8toi4 (yhead(1), iheadv(1),8)
CALL pbread(kfilet,yhead(1),8*8,kret)
IF (kret/=8*8) CALL finish('ReadOneBlock','nudging data inconsistent TEM')
CALL util_i8toi4 (yhead(1), iheadt(1),8)
975:
! check header consistency
IF ((iheadd(3)==iheadv(3)).AND.(iheadd(4)==iheadv(4)) &
& .AND.(iheadd(3)==iheadt(3)).AND.(iheadd(4)==iheadt(4))) THEN
!
980: ! calculate corresponding time step
iday = iymd2c(iheadd(3))
ndgstep2 = ((iday-ncbase)*dayl+iheadd(4)*3600-ntbase)/dtime
ihead2d = ic2ymd(iday) ! YYMMDD
985: ihead2t = iheadd(4)*10000 ! HH0000 (hhmmss)
!
! read data block
#ifdef CRAY
990: CALL pbread(kfiled,phbuf,ilens1*8,kret)
IF (kret/=ilens1*8) CALL finish('ReadOneBlock','nudging data fault DIV')
iamount = ilens1
CALL util_reshape(sdobs2,phbuf,iamount)
995: CALL pbread(kfilev,phbuf,ilens1*8,kret)
IF (kret/=ilens1*8) CALL finish('ReadOneBlock','nudging data fault VOR')
iamount = ilens1
CALL util_reshape(svobs2,phbuf,iamount)
1000: CALL pbread(kfilet,phbuf,ilens2*8,kret)
IF (kret/=ilens2*8) CALL finish('ReadOneBlock','nudging data fault TEM')
iamount = ilens2
CALL util_reshape(stobs2,phbuf,iamount)
#else
1005: CALL pbread(kfiled,phbuf,ilens1*8,kret)
IF (kret/=ilens1*8) CALL finish('ReadOneBlock','nudging data fault DIV')
CALL util_cray2ieee(phbuf,zhbuf,ilens1)
sdobs2 = RESHAPE(zhbuf,(/nlev,2,nsp/))
1010: CALL pbread(kfilev,phbuf,ilens1*8,kret)
IF (kret/=ilens1*8) CALL finish('ReadOneBlock','nudging data fault VOR')
CALL util_cray2ieee(phbuf,zhbuf,ilens1)
svobs2 = RESHAPE(zhbuf,(/nlev,2,nsp/))
1015: CALL pbread(kfilet,phbuf,ilens2*8,kret)
IF (kret/=ilens2*8) CALL finish('ReadOneBlock','nudging data fault TEM')
CALL util_cray2ieee(phbuf,zhbuf,ilens2)
stobs2 = RESHAPE(zhbuf,(/nlevp1,2,nsp/))
#endif
1020:
IF (lnudgdbx) WRITE(nout,*) &
& 'ReadOneBlock: at date= ',ihead2d &
& ,' time= ',ihead2t &
& ,' refstep= ',ndgstep2
1025:
!! if (lnmi) call NMI_Make(1)
EXIT ! read block OBS2 HEAD2 success
ELSE
CALL finish('ReadOneBlock','header synchronisation fault')
1030: ENDIF
ENDDO
#ifndef CRAY
1035: DEALLOCATE (zhbuf)
#endif
DEALLOCATE (phbuf)
RETURN
1040:
END SUBROUTINE ReadOneBlock
!======================================================================
SUBROUTINE OpenOneBlock
1045:
! open nudging data block
INTEGER :: kret
IF ((0 < ndgblock).AND.(ndgblock < 4)) THEN
1050:
! open nudging data block
WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(3*(ndgblock-1)+1)
CALL pbopen(kfiled,cfile,'r',kret)
1055: IF(kret/=0) CALL finish('OpenOneBlock','data file available? DIV')
WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(3*(ndgblock-1)+2)
CALL pbopen(kfilev,cfile,'r',kret)
IF(kret/=0) CALL finish('OpenOneBlock','data file available? VOR')
1060:
WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(3*(ndgblock-1)+3)
CALL pbopen(kfilet,cfile,'r',kret)
IF(kret/=0) CALL finish('OpenOneBlock','data file available? TEM')
1065: #ifndef CRAY
WRITE(nout,*) ' Attention unit ', &
& ndunit(3*(ndgblock-1)+1),&
& ndunit(3*(ndgblock-1)+2),&
& ndunit(3*(ndgblock-1)+3),' : convert NUDGE data to IEEE'
1070: #endif
ELSE
CALL finish('OpenOneBlock','nudging data block number wrong')
ENDIF
1075:
RETURN
END SUBROUTINE OpenOneBlock
!======================================================================
1080:
SUBROUTINE CloseBlock
! close nudging data block
INTEGER :: kret
CALL pbclose(kfiled,kret)
1085: CALL pbclose(kfilev,kret)
CALL pbclose(kfilet,kret)
END SUBROUTINE CloseBlock
!======================================================================
1090: SUBROUTINE NudgingReadSST
! Read ECHAM SST for nudging experiments
! call every time step in time step loop
!
! Local scalars
1095: INTEGER :: krtim, kdag, iday, id, im, iy, ihh, kret, jr, ipos
REAL :: zsec
! Local arrays
REAL, ALLOCATABLE :: zhbuf(:)
1100: INTEGER :: ihead(8)
CHARACTER (8) :: yhead(8)
! Intrinsic functions
INTRINSIC MOD
1105:
! Executable statements
IF (ltwin) RETURN
1110: ! read next SST set
IF (nsstinc/=0) THEN
! ---------------------------------------------------------------
1115: ! KHOURS gives interval where SSTs are available.
! KHOURS=24 -> SST available once a day,
! KHOURS=12 -> SST available twice a day (only for twin exp.)
! NSSTINC setup in RUNCTL
! ---------------------------------------------------------------
1120:
krtim = ntbase + nstep*dtime - nsstoff*3600 ! sst time step offset in sec
kdag = nsstinc*3600 ! sst time increment in sec
lsstn = (MOD(krtim,kdag) == 0).OR.(iheads == 0)
1125: IF(lsstn) THEN
! new sst field neccessary
!
! Read one set of SSTs
! Calculate the model time
1130: IF (iheads==0) THEN
! get last sst date/time befor actual time step
zsec = ntbase+dtime*nstep - MOD(krtim,kdag)
ELSE
zsec = ntbase+dtime*nstep
1135: ENDIF
!
! calculate excpeted date and time
iday = ncbase + zsec/dayl + 0.0001
CALL cd2dat(iday,id,im,iy)
1140: ihh = INT(MOD(zsec,dayl)/3600)
iheads = 1000000*iy+10000*im+100*id+ihh
iheads = MOD(iheads,100000000) ! remove century
IF (lnudgdbx) WRITE(nout,*) &
1145: & 'NudgingReadSST: need DATE/TIME ??YYMMDDhh ',iheads
IF (sstblock<0) THEN
! set and open start block
sstblock = -sstblock
1150: WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(9+sstblock)
CALL pbopen(kfiles,cfile,'r',kret)
IF (kret/=0) CALL finish('NudgingReadSST','nudging sst files available?')
ipos = 0
ENDIF
1155:
! read next sst record
DO
! read header
CALL pbread(kfiles,yhead(1),8*8,kret)
1160: IF (kret/=8*8) THEN
CALL pbclose(kfiles,kret)
sstblock = sstblock+1
IF (sstblock>3) CALL finish ('NudgingReadSST','need more sst files')
WRITE(cfile,'(a5,i2.2)') 'fort.',ndunit(9+sstblock)
1165: CALL pbopen(kfiles,cfile,'r',kret)
WRITE (nout,*) 'NudgingReadSST open next SST data file ',sstblock,ndunit(9+sstblock)
ipos = 0
CYCLE
ENDIF
1170: CALL util_i8toi4(yhead(1), ihead(1), 8)
! check header without century
ihead(3) = MOD(ihead(3),1000000)
IF (ihead(3)*100+ihead(4) > iheads) CALL finish('NudgingReadSST','sst date fault')
IF (ihead(3)*100+ihead(4) == iheads) EXIT
1175: ! skip data set
ipos = ipos + 8 + nlon*ngl
CALL pbseek(kfiles,ipos*8,0,kret)
WRITE (nout,*) 'NudgingReadSST: skip SST ... header ',ihead
ENDDO
1180: ! found record
#ifndef CRAY
WRITE(*,*) ' Attention convert SST data to IEEE'
ALLOCATE(zhbuf(nlon))
1185: zhbuf(:) = 0.0
#endif
! read odd latitudes
DO jr=1,ngl,2
1190: CALL pbread(kfiles,sstn(1,jr),nlon*8,kret)
#ifndef CRAY
CALL util_cray2ieee(sstn(1,jr),zhbuf,nlon)
sstn(:,jr)=zhbuf(:)
#endif
1195: ENDDO
! read even latitudes
DO jr=ngl,2,-2
CALL pbread(kfiles,sstn(1,jr),nlon*8,kret)
#ifndef CRAY
1200: CALL util_cray2ieee(sstn(1,jr),zhbuf,nlon)
sstn(:,jr)=zhbuf(:)
#endif
ENDDO
1205: #ifndef CRAY
DEALLOCATE (zhbuf)
#endif
IF (lnudgdbx) WRITE(nout,*) &
1210: & 'NudgingReadSST: read new SST ... header ',ihead
ENDIF
ENDIF
1215: RETURN
END SUBROUTINE NudgingReadSST
!======================================================================
1220: SUBROUTINE NudgingSSTnew
!
! call every time step in latitude loop
!
!-------------------------------------------------------
1225: ! Read NMC-sea surface temperature analysis retrieved
! from the mars archive, instead of the climatology provided
! for ECHAM. SST's are available every 24 hours
!--------------------------------------------------------
!
1230: ! Local scalars
INTEGER :: irow, jrow, jn
REAL, PARAMETER :: zdt=0.01 ! correct temperature near freezing point
REAL :: zts
1235: ! Intrinsic functions
INTRINSIC MAX, MIN
! Executable statements
1240: IF (nsstinc==0.OR.ltwin) THEN
! use standard echam sst
CALL clsst
1245: IF (ltwin) sstn(:,:) = tsm(:,:)
ELSE
! update sea surface temperatures at every time step
1250:
irow = nrow(1)
jrow = nrow(2)
DO jn = 1,nlon
1255:
zts = sstn(jn,irow)
IF (slmm(jn,jrow)<=0.5) THEN
! sea point
1260: IF (zts<=ctfreez) THEN
! below freezing level
tsm(jn,jrow) =MIN(zts,ctfreez-zdt)
tsm1m(jn,jrow)=MIN(zts,ctfreez-zdt)
ELSE
1265: ! above freezing level
tsm(jn,jrow) =MAX(zts,ctfreez+zdt)
tsm1m(jn,jrow)=MAX(zts,ctfreez+zdt)
ENDIF
ENDIF
1270:
ENDDO
IF (lsstn.AND..NOT.ltwin) THEN
! read new sst field at this timestep
1275: ! correction of seaice skin-temperature
auxil1m(:,jrow)=tsm(:,jrow)
auxil2m(:,jrow)=tsm(:,jrow)
IF (lnudgdbx.AND.(irow==ngl)) WRITE(nout,*) &
1280: & 'NudgingSSTnew: SEAICE correction at NSTEP= ',nstep
END IF
ENDIF
1285:
RETURN
END SUBROUTINE NudgingSSTnew
!======================================================================
1290:
SUBROUTINE NudgingOut
! Write out instantaneous (OUTTINS)
! and accumulated (OUTTACC) tendencies
1295: REAL :: tfact
INTEGER :: jlev, ierr, iret
REAL :: worksp(2,nsp)
1300: IF (ltwin) RETURN
! define GRIB block 1
ksec1(1) = nudging_table
1305:
CALL set_output_time
ksec1(10) = year
ksec1(11) = month
1310: ksec1(12) = day
ksec1(13) = hour
ksec1(14) = minute
ksec1(21) = century
1315: level_type = 109
ksec1(7) = level_type
ksec4(1) = n2sp
1320: WRITE (nout,'(a,i2.2,a1,i2.2,2x,i2.2,a1,i2.2,a1,i2.2,i2.2)') &
& 'NudgingOut: tendencies output at ... ', &
& hour, ':', minute, day, '.', month, '.', century-1, year
! store only nudging part
1325:
! 1. write instantaneous values
IF (lppp) THEN ! 1.1 Pressure
code_parameter = ndgcode(1)
1330: ksec1(6) = code_parameter
level_p2 = 0
ksec1(9) = level_p2
worksp(:,:) = flagn(nlevp1,:,:)*stten(nlevp1,:,:)
CALL NudgStoreSP
1335: ENDIF
IF (lppt) THEN ! 1.2 Temperature
code_parameter = ndgcode(2)
ksec1(6) = code_parameter
1340: DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = flagn(jlev,:,:)*stten(jlev,:,:)
CALL NudgStoreSP
1345: END DO
END IF
IF (lppd) THEN ! 1.3 Divergence
code_parameter = ndgcode(3)
1350: ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = flagn(jlev,:,:)*sdten(jlev,:,:)
1355: CALL NudgStoreSP
END DO
END IF
IF (lppvo) THEN ! 1.4 Vorticity
1360: code_parameter = ndgcode(4)
ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
1365: worksp(:,:) = flagn(jlev,:,:)*svten(jlev,:,:)
CALL NudgStoreSP
END DO
END IF
1370: ! 2. accumulated data
! proportionality factor, output given in VALUE/sec
IF ( (lnudgini.AND.lres.AND.nstep==1) .OR. &
& (lnudgini.AND.(.NOT.lres).AND.nstep==0) ) THEN
1375: ! neccessary for initialization
tfact = 1.0
ELSE
tfact = 1.0/nptime
ENDIF
1380:
IF (lnudgdbx) WRITE(nout,*) &
& ' correct accumulated data with = ',tfact
IF (lppp) THEN ! 2.1 Pressure
1385: code_parameter = ndgcode(5)
ksec1(6) = code_parameter
level_p2 = 0
ksec1(9) = level_p2
worksp(:,:) = flagn(nlevp1,:,:)*sttac(nlevp1,:,:)*tfact
1390: CALL NudgStoreSP
ENDIF
IF (lppt) THEN ! 2.2 Temperature
code_parameter = ndgcode(6)
1395: ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = flagn(jlev,:,:)*sttac(jlev,:,:)*tfact
1400: CALL NudgStoreSP
END DO
END IF
IF (lppd) THEN ! 2.3 Divergence
1405: code_parameter = ndgcode(7)
ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
1410: worksp(:,:) = flagn(jlev,:,:)*sdtac(jlev,:,:)*tfact
CALL NudgStoreSP
END DO
END IF
1415: IF (lppvo) THEN ! 2.4 Vorticity
code_parameter = ndgcode(8)
ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
1420: ksec1(9) = level_p2
worksp(:,:) = flagn(jlev,:,:)*svtac(jlev,:,:)*tfact
CALL NudgStoreSP
END DO
END IF
1425:
IF (lnudgwobs) THEN
! store diagnostic part
1430:
! 3 write instantaneous values
IF (lppp) THEN ! 3.1 Pressure
code_parameter = ndgcode(9)
1435: ksec1(6) = code_parameter
level_p2 = 0
ksec1(9) = level_p2
worksp(:,:) = flago(nlevp1,:,:)*stten(nlevp1,:,:)
CALL NudgStoreSP
1440: ENDIF
IF (lppt) THEN ! 3.2 Temperature
code_parameter = ndgcode(10)
ksec1(6) = code_parameter
1445: DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = flago(jlev,:,:)*stten(jlev,:,:)
CALL NudgStoreSP
1450: END DO
END IF
IF (lppd) THEN ! 3.3 Divergence
code_parameter = ndgcode(11)
1455: ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = flago(jlev,:,:)*sdten(jlev,:,:)
1460: CALL NudgStoreSP
END DO
END IF
IF (lppvo) THEN ! 3.4 Vorticity
1465: code_parameter = ndgcode(12)
ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
1470: worksp(:,:) = flago(jlev,:,:)*svten(jlev,:,:)
CALL NudgStoreSP
END DO
END IF
1475: ! 4. accumulated data
IF (lppp) THEN ! 4.1 Pressure
code_parameter = ndgcode(13)
ksec1(6) = code_parameter
1480: level_p2 = 0
ksec1(9) = level_p2
worksp(:,:) = flago(nlevp1,:,:)*sttac(nlevp1,:,:)*tfact
CALL NudgStoreSP
ENDIF
1485:
IF (lppt) THEN ! 4.2 Temperature
code_parameter = ndgcode(14)
ksec1(6) = code_parameter
DO jlev = 1, nlev
1490: level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = flago(jlev,:,:)*sttac(jlev,:,:)*tfact
CALL NudgStoreSP
END DO
1495: END IF
IF (lppd) THEN ! 4.3 Divergence
code_parameter = ndgcode(15)
ksec1(6) = code_parameter
1500: DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = flago(jlev,:,:)*sdtac(jlev,:,:)*tfact
CALL NudgStoreSP
1505: END DO
END IF
IF (lppvo) THEN ! 4.4 Vorticity
code_parameter = ndgcode(16)
1510: ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = flago(jlev,:,:)*svtac(jlev,:,:)*tfact
1515: CALL NudgStoreSP
END DO
END IF
! 5 write observations only instantaneous values
1520:
IF (lppp) THEN ! 5.1 Pressure
code_parameter = ndgcode(17)
ksec1(6) = code_parameter
level_p2 = 0
1525: ksec1(9) = level_p2
worksp(:,:) = stobs(nlevp1,:,:)
CALL NudgStoreSP
ENDIF
1530: IF (lppt) THEN ! 5.2 Temperature
code_parameter = ndgcode(18)
ksec1(6) = code_parameter
DO jlev = 1, nlev
level_p2 = jlev
1535: ksec1(9) = level_p2
worksp(:,:) = stobs(jlev,:,:)
CALL NudgStoreSP
END DO
END IF
1540:
IF (lppd) THEN ! 3.3 Divergence
code_parameter = ndgcode(19)
ksec1(6) = code_parameter
DO jlev = 1, nlev
1545: level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = sdobs(jlev,:,:)
CALL NudgStoreSP
END DO
1550: END IF
IF (lppvo) THEN ! 3.4 Vorticity
code_parameter = ndgcode(20)
ksec1(6) = code_parameter
1555: DO jlev = 1, nlev
level_p2 = jlev
ksec1(9) = level_p2
worksp(:,:) = svobs(jlev,:,:)
CALL NudgStoreSP
1560: END DO
END IF
END IF
1565: sttac(:,:,:) = 0.0
sdtac(:,:,:) = 0.0
svtac(:,:,:) = 0.0
RETURN
1570:
CONTAINS
SUBROUTINE NudgStoreSP
1575: ! External subroutines
EXTERNAL codegb5, pbwrite
#ifdef EMOS
CALL gribex (ksec0, ksec1, ksec2_sp, psec2, ksec3, psec3, &
1580: & ksec4, worksp, n2sp, kgrib, kleng, kword, 'C', ierr)
#else
CALL codegb5 (worksp,n2sp,16,nbit,ksec1,ksec2_sp,vct,2*nlevp1, &
& kgrib,klengo,kword,0,ierr)
#endif
1585: CALL pbwrite(nunitdf,kgrib,kword*iobyte,iret)
IF (iret /= kword*iobyte) &
& CALL finish('NudgeStoreSP', &
& 'I/O error on spectral harmonics output - disk full?')
1590: RETURN
END SUBROUTINE NudgStoreSP
END SUBROUTINE NudgingOut
!======================================================================
1595:
SUBROUTINE WriteTwin
INTEGER :: isecs, ihead(8), ilens1, ilens2, iilen, kret, jr
CHARACTER (8) :: yhead(8)
REAL, ALLOCATABLE :: zhbuf(:)
1600:
! get time at the day
isecs = ntbase + (nstep+1)*dtime
! WRITE (nout,*) 'WriteTwin:: secs at the day ',isecs
! check output time
1605:
IF (MOD(isecs-ntwoff*3600,ntwinc*3600)==0) THEN
CALL set_output_time
1610:
ilens1 = n2sp*nlev
ilens2 = n2sp*nlevp1
iilen = MAX(ilens1,ilens2)
1615: ALLOCATE(zhbuf(iilen))
zhbuf(:) = 0.0
WRITE(nout,*) 'WriteTwin:: store twin data at ',isecs
1620: ihead(2) = nlev
ihead(3) = (century-1)*1000000+year*10000+month*100+day ! YYMMDD
ihead(4) = hour ! HH
ihead(5) = n2sp
ihead(6) = nlev
1625: ihead(7) = 0
ihead(8) = 0
! store vorticity
ihead(1) = 138
1630: CALL util_i4toi8(ihead(1),yhead(1),8)
CALL pbwrite(kfilev,yhead(1),8*8,kret)
#ifdef CRAY
CALL util_reshape(zhbuf,svo,ilens1)
CALL pbwrite(kfilev,zhbuf,ilens1*8,kret)
1635: #else
CALL util_ieee2cray(svo,zhbuf,ilens1)
CALL pbwrite(kfilev,zhbuf,ilens1*8,kret)
#endif
1640: ! store divergence
ihead(1) = 155
CALL util_i4toi8(ihead(1),yhead(1),8)
CALL pbwrite(kfiled,yhead(1),8*8,kret)
#ifdef CRAY
1645: CALL util_reshape(zhbuf,sd,ilens1)
CALL pbwrite(kfiled,zhbuf,ilens1*8,kret)
#else
CALL util_ieee2cray(sd,zhbuf,ilens1)
CALL pbwrite(kfiled,zhbuf,ilens1*8,kret)
1650: #endif
! store temperature and pressure
ihead(1) = 130
ihead(6) = nlevp1
1655: CALL util_i4toi8(ihead(1),yhead(1),8)
CALL pbwrite(kfilet,yhead(1),8*8,kret)
#ifdef CRAY
CALL util_reshape(zhbuf,stp,ilens2)
CALL pbwrite(kfilet,zhbuf,ilens2*8,kret)
1660: #else
CALL util_ieee2cray(stp,zhbuf,ilens2)
CALL pbwrite(kfilet,zhbuf,ilens2*8,kret)
#endif
1665: ! store SST
ihead(1) = 169
CALL util_i4toi8(ihead(1),yhead(1),8)
CALL pbwrite(kfiles,yhead(1),8*8,kret)
! write odd latitudes
1670: DO jr=1,ngl,2
#ifndef CRAY
CALL util_ieee2cray(zhbuf,sstn(1,jr),nlon)
#else
zhbuf(1:nlon)= sstn(:,jr)
1675: #endif
CALL pbwrite(kfiles,zhbuf(1),nlon*8,kret)
ENDDO
! write even latitudes
1680: DO jr=ngl,2,-2
#ifndef CRAY
CALL util_ieee2cray(zhbuf,sstn(1,jr),nlon)
#else
zhbuf(1:nlon)=sstn(:,jr)
1685: #endif
CALL pbwrite(kfiles,zhbuf(1),nlon*8,kret)
ENDDO
ENDIF
1690:
RETURN
END SUBROUTINE WriteTwin
!======================================================================
1695:
SUBROUTINE NudgingRerun(flag)
! store and restore accumulated fields for nudging
INTEGER :: flag, rrunit, iistep, iinsp, iinlev
1700:
IF (ltwin) RETURN
rrunit = ndunit(13)
1705: WRITE(nout,*) 'NudgingRerun: unit ',rrunit
SELECT CASE(flag)
CASE(-1) ! store restart data
1710: REWIND(rrunit)
WRITE(rrunit) nstep+1,nsp,nlev
WRITE(rrunit) sdtac
WRITE(rrunit) svtac
WRITE(rrunit) sttac
1715: WRITE(nout,*) ' store accu-tend at NSTEP+1 = ',nstep+1
CLOSE(rrunit)
CASE(1) ! reload restart data
1720: sdtac(:,:,:) = 0.0
svtac(:,:,:) = 0.0
sttac(:,:,:) = 0.0
REWIND(rrunit)
READ(rrunit,err=100) iistep,iinsp,iinlev
1725: IF ( (iinsp/=nsp).OR.(iinlev/=nlev) ) THEN
WRITE(nout,*) 'NudgingRerun: ERROR dimension in nudging data rerunfile'
GOTO 100
END IF
READ(rrunit) sdtac
1730: READ(rrunit) svtac
READ(rrunit) sttac
WRITE(nout,*) 'NudgingRerun: read data STEP ',iistep,' model step is ',nstep
CLOSE(rrunit)
1735: CASE default
WRITE(nout,*) ' set accu-tend to ZERO'
sdtac(:,:,:) = 0.0
svtac(:,:,:) = 0.0
sttac(:,:,:) = 0.0
1740: END SELECT
100 CONTINUE
RETURN
1745:
END SUBROUTINE NudgingRerun
END MODULE mo_nudging
Info Section
uses: all, block, explicit, implicit, mo_constants
mo_control, mo_doctor, mo_exception, mo_grib, mo_kind
mo_memory_g3a, mo_memory_sp, mo_mpi, mo_nmi, mo_nudging_buffer
mo_physc2, mo_post, mo_start_dataset, mo_truncation, mo_year
standard, time, truncation
includes: ndgctl.inc
calls: cd2dat, clockm, closeblock, clsst, codegb5
finish, getnudgedata, gribex, nmi_make, nudgingrerun
nudgstoresp, openoneblock, pbclose, pbopen, pbread
pbseek, pbwrite, readoneblock, set_output_time, util_cray2ieee
util_i4toi8, util_i8toi4, util_ieee2cray, util_reshape, writetwin
back to top
ECHAM 4 vf90 (C) 1998 Max-Planck-Institut für Meteorologie, Hamburg
Wed Nov 24 01:25:21 CST 1999
HTML derived from FORTRAN source by f2html.pl v0.3 (C) 1997,98 Beroud Jean-Marc.