global_cycle  1.13.0
cycle.f90
Go to the documentation of this file.
1 
4 
111  PROGRAM sfc_drv
112 
113  use mpi
114 
115  IMPLICIT NONE
116 !
117  CHARACTER(LEN=3) :: donst
118  INTEGER :: idim, jdim, lsoil, lugb, iy, im, id, ih, ialb
119  INTEGER :: isot, ivegsrc, lensfc, zsea1_mm, zsea2_mm, ierr
120  INTEGER :: nprocs, myrank, num_threads, num_parthds, max_tasks
121  REAL :: fh, deltsfc, zsea1, zsea2
122  LOGICAL :: use_ufo, do_nsst, do_lndinc, do_sfccycle, frac_grid
123 !
124  NAMELIST/namcyc/ idim,jdim,lsoil,lugb,iy,im,id,ih,fh,&
125  deltsfc,ialb,use_ufo,donst, &
126  do_sfccycle,isot,ivegsrc,zsea1_mm, &
127  zsea2_mm, max_tasks, do_lndinc, frac_grid
128 !
129  DATA idim,jdim,lsoil/96,96,4/
130  DATA iy,im,id,ih,fh/1997,8,2,0,0./
131  DATA lugb/51/, deltsfc/0.0/, ialb/1/, max_tasks/99999/
132  DATA isot/1/, ivegsrc/2/, zsea1_mm/0/, zsea2_mm/0/
133 !
134  CALL mpi_init(ierr)
135  CALL mpi_comm_size(mpi_comm_world, nprocs, ierr)
136  CALL mpi_comm_rank(mpi_comm_world, myrank, ierr)
137 
138  if (myrank==0) call w3tagb('GLOBAL_CYCLE',2018,0179,0055,'NP20')
139 
140  num_threads = num_parthds()
141 
142  print*
143  print*,"STARTING CYCLE PROGRAM ON RANK ", myrank
144  print*,"RUNNING WITH ", nprocs, "TASKS"
145  print*,"AND WITH ", num_threads, " THREADS."
146 
147  use_ufo = .false.
148  donst = "NO"
149  do_lndinc = .false.
150  do_sfccycle = .true.
151  frac_grid = .false.
152 
153  print*
154  print*,"READ NAMCYC NAMELIST."
155 
156  CALL baopenr(36, "fort.36", ierr)
157  IF (ierr /= 0) THEN
158  print*,'FATAL ERROR READING FORT.36 NAMELIST. IERR: ', ierr
159  CALL mpi_abort(mpi_comm_world, 32, ierr)
160  ENDIF
161  READ(36, nml=namcyc, iostat=ierr)
162  IF (ierr /= 0) THEN
163  print*,'FATAL ERROR READING FORT.36 NAMELIST. IERR: ', ierr
164  CALL mpi_abort(mpi_comm_world, 33, ierr)
165  ENDIF
166 !IF (MYRANK==0) WRITE(6,NAMCYC)
167 
168  IF (max_tasks < 99999 .AND. myrank > (max_tasks - 1)) THEN
169  print*,"USER SPECIFIED MAX NUMBER OF TASKS: ", max_tasks
170  print*,"WILL NOT RUN CYCLE PROGRAM ON RANK: ", myrank
171  GOTO 333
172  ENDIF
173 
174  lensfc = idim*jdim ! TOTAL NUMBER OF POINTS FOR THE CUBED-SPHERE TILE
175 
176  zsea1 = float(zsea1_mm) / 1000.0 ! CONVERT FROM MM TO METERS
177  zsea2 = float(zsea2_mm) / 1000.0
178 
179  IF (donst == "YES") THEN
180  do_nsst=.true.
181  ELSE
182  do_nsst=.false.
183  ENDIF
184 
185  print*
186  IF (myrank==0) print*,"LUGB,IDIM,JDIM,ISOT,IVEGSRC,LSOIL,DELTSFC,IY,IM,ID,IH,FH: ", &
187  lugb,idim,jdim,isot,ivegsrc,lsoil,deltsfc,iy,im,id,ih,fh
188 
189  CALL sfcdrv(lugb,idim,jdim,lensfc,lsoil,deltsfc, &
190  iy,im,id,ih,fh,ialb, &
191  use_ufo,do_nsst,do_sfccycle,do_lndinc, &
192  frac_grid,zsea1,zsea2,isot,ivegsrc,myrank)
193 
194  print*
195  print*,'CYCLE PROGRAM COMPLETED NORMALLY ON RANK: ', myrank
196 
197  333 CONTINUE
198 
199  CALL mpi_barrier(mpi_comm_world, ierr)
200 
201  if (myrank==0) call w3tage('GLOBAL_CYCLE')
202 
203  CALL mpi_finalize(ierr)
204 
205  stop
206 
207  END PROGRAM sfc_drv
208 
317  SUBROUTINE sfcdrv(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, &
318  IY,IM,ID,IH,FH,IALB, &
319  USE_UFO,DO_NSST,DO_SFCCYCLE,DO_LNDINC,&
320  FRAC_GRID,ZSEA1,ZSEA2,ISOT,IVEGSRC,MYRANK)
321 !
322  USE read_write_data
323  use machine
324  USE mpi
325  USE land_increments, ONLY: gaussian_to_fv3_interp, &
326  add_increment_soil, &
327  add_jedi_increment_snow, &
328  calculate_landinc_mask, &
329  apply_land_da_adjustments_soil, &
330  apply_land_da_adjustments_snd, &
331  lsm_noah, lsm_noahmp
332 
333  IMPLICIT NONE
334 
335  INTEGER, INTENT(IN) :: IDIM, JDIM, LENSFC, LSOIL, IALB
336  INTEGER, INTENT(IN) :: LUGB, IY, IM, ID, IH
337  INTEGER, INTENT(IN) :: ISOT, IVEGSRC, MYRANK
338 
339  LOGICAL, INTENT(IN) :: USE_UFO, DO_NSST,DO_SFCCYCLE
340  LOGICAL, INTENT(IN) :: DO_LNDINC, FRAC_GRID
341 
342  REAL, INTENT(IN) :: FH, DELTSFC, ZSEA1, ZSEA2
343 
344  INTEGER, PARAMETER :: NLUNIT=35
345  INTEGER, PARAMETER :: SZ_NML=1
346 
347  CHARACTER(LEN=5) :: TILE_NUM
348  CHARACTER(LEN=500) :: NST_FILE
349  CHARACTER(LEN=500) :: GSI_SOI_FILE,JEDI_SOI_FILE,JEDI_SNO_FILE
350  CHARACTER(LEN=4) :: INPUT_NML_FILE(SZ_NML)
351 
352  INTEGER :: I, IERR
353  INTEGER :: I_INDEX(LENSFC), J_INDEX(LENSFC)
354  INTEGER :: IDUM(IDIM,JDIM)
355  integer :: num_parthds, num_threads
356 
357  LOGICAL :: IS_NOAHMP
358  INTEGER :: LSM
359 
360  real(kind=kind_io8) :: min_ice(lensfc)
361 
362  REAL :: SLMASK(LENSFC), OROG(LENSFC)
363  REAL :: SIHFCS(LENSFC), SICFCS(LENSFC)
364  REAL :: SITFCS(LENSFC), TSFFCS(LENSFC)
365  REAL :: SWEFCS(LENSFC), ZORFCS(LENSFC)
366  REAL :: ALBFCS(LENSFC,4), TG3FCS(LENSFC)
367  REAL :: CNPFCS(LENSFC), SMCFCS(LENSFC,LSOIL)
368  REAL :: STCFCS(LENSFC,LSOIL), SLIFCS(LENSFC)
369  REAL :: AISFCS(LENSFC), F10M(LENSFC)
370  REAL :: VEGFCS(LENSFC), VETFCS(LENSFC)
371  REAL :: SOTFCS(LENSFC), ALFFCS(LENSFC,2)
372  REAL :: CVFCS(LENSFC), CVTFCS(LENSFC)
373  REAL :: CVBFCS(LENSFC), TPRCP(LENSFC)
374  REAL :: SRFLAG(LENSFC), SNDFCS(LENSFC)
375  REAL :: SLCFCS(LENSFC,LSOIL), VMXFCS(LENSFC)
376  REAL :: VMNFCS(LENSFC), T2M(LENSFC)
377  REAL :: Q2M(LENSFC), SLPFCS(LENSFC)
378  REAL :: ABSFCS(LENSFC), OROG_UF(LENSFC)
379  REAL :: USTAR(LENSFC), SOCFCS(LENSFC)
380  REAL :: FMM(LENSFC), FHH(LENSFC)
381  REAL :: RLA(LENSFC), RLO(LENSFC)
382  REAL(KIND=4) :: ZSOIL(LSOIL)
383  REAL :: SIG1T(LENSFC)
386  REAL, ALLOCATABLE :: STC_BCK(:,:), SMC_BCK(:,:), SLC_BCK(:,:)
387  REAL, ALLOCATABLE :: SLIFCS_FG(:), SICFCS_FG(:)
388  INTEGER, ALLOCATABLE :: LANDINC_MASK_FG(:), LANDINC_MASK(:)
389  REAL, ALLOCATABLE :: SND_BCK(:), SND_INC(:), SWE_BCK(:)
390  REAL(KIND=KIND_IO8), ALLOCATABLE :: SLMASKL(:), SLMASKW(:), LANDFRAC(:)
391 
392  TYPE(NSST_DATA) :: NSST
393  real, dimension(idim,jdim) :: tf_clm,tf_trd,sal_clm
394  real, dimension(lensfc) :: tf_clm_tile,tf_trd_tile,sal_clm_tile
395  INTEGER :: veg_type_landice
396  INTEGER, DIMENSION(LENSFC) :: STC_UPDATED, SLC_UPDATED
397  REAL, DIMENSION(LENSFC,LSOIL) :: STCINC, SLCINC
398 
399  LOGICAL :: FILE_EXISTS, DO_SOI_INC_GSI, DO_SOI_INC_JEDI, DO_SNO_INC_JEDI
400  CHARACTER(LEN=3) :: RANKCH
401  INTEGER :: lsoil_incr
402 
403 !--------------------------------------------------------------------------------
404 ! NST_FILE is the path/name of the gaussian GSI file which contains NSST
405 ! increments.
406 !--------------------------------------------------------------------------------
407 
408  NAMELIST/namsfcd/ nst_file, lsoil_incr, do_sno_inc_jedi, do_soi_inc_jedi, do_soi_inc_gsi
409 
410  DATA nst_file/'NULL'/
411 
412  do_sno_inc_jedi = .false.
413  do_soi_inc_gsi = .false.
414  do_soi_inc_jedi = .false.
415  lsoil_incr = 3 !default
416 
417 
418  sig1t = 0.0 ! Not a dead start!
419 
420  input_nml_file = "NULL"
421 
422  CALL baopenr(37, "fort.37", ierr)
423  IF (ierr /= 0) THEN
424  print*,'FATAL ERROR OPENING FORT.37 NAMELIST. IERR: ', ierr
425  CALL mpi_abort(mpi_comm_world, 30, ierr)
426  ENDIF
427  READ (37, nml=namsfcd, iostat=ierr)
428  IF (ierr /= 0) THEN
429  print*,'FATAL ERROR READING FORT.37 NAMELIST. IERR: ', ierr
430  CALL mpi_abort(mpi_comm_world, 31, ierr)
431  ENDIF
432 
433  print*
434  print*,'IN ROUTINE SFCDRV,IDIM=',idim,'JDIM=',jdim,'FH=',fh
435 
436 !--------------------------------------------------------------------------------
437 ! READ THE OROGRAPHY AND GRID POINT LAT/LONS FOR THE CUBED-SPHERE TILE.
438 !--------------------------------------------------------------------------------
439 
440  ALLOCATE(landfrac(lensfc))
441  IF(frac_grid) THEN
442  print*,'- RUNNING WITH FRACTIONAL GRID.'
443  CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc,landfrac=landfrac)
444  ELSE
445  CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc)
446  landfrac=-999.9
447  ENDIF
448 
449  DO i = 1, idim
450  idum(i,:) = i
451  ENDDO
452 
453  i_index = reshape(idum, (/lensfc/))
454 
455  DO i = 1, jdim
456  idum(:,i) = i
457  ENDDO
458 
459  j_index = reshape(idum, (/lensfc/) )
460 
461  IF (do_nsst) THEN
462  print*
463  print*,"WILL PROCESS NSST RECORDS."
464  ALLOCATE(nsst%C_0(lensfc))
465  ALLOCATE(nsst%C_D(lensfc))
466  ALLOCATE(nsst%D_CONV(lensfc))
467  ALLOCATE(nsst%DT_COOL(lensfc))
468  ALLOCATE(nsst%IFD(lensfc))
469  ALLOCATE(nsst%QRAIN(lensfc))
470  ALLOCATE(nsst%TREF(lensfc))
471  ALLOCATE(nsst%TFINC(lensfc))
472  ALLOCATE(nsst%W_0(lensfc))
473  ALLOCATE(nsst%W_D(lensfc))
474  ALLOCATE(nsst%XS(lensfc))
475  ALLOCATE(nsst%XT(lensfc))
476  ALLOCATE(nsst%XTTS(lensfc))
477  ALLOCATE(nsst%XU(lensfc))
478  ALLOCATE(nsst%XV(lensfc))
479  ALLOCATE(nsst%XZ(lensfc))
480  ALLOCATE(nsst%XZTS(lensfc))
481  ALLOCATE(nsst%Z_C(lensfc))
482  ALLOCATE(nsst%ZM(lensfc))
483  ALLOCATE(slifcs_fg(lensfc))
484  ALLOCATE(sicfcs_fg(lensfc))
485  ENDIF
486 
487 IF (do_lndinc) THEN
488  ! identify variables to be updated, and allocate arrays.
489  IF (do_soi_inc_gsi .and. do_soi_inc_jedi) THEN
490  print*
491  print*, 'FATAL ERROR: Can not do gsi and jedi soil updates at the same time, choose one!'
492  CALL mpi_abort(mpi_comm_world, 15, ierr)
493  ENDIF
494  IF (do_soi_inc_gsi .or. do_soi_inc_jedi) THEN
495  print*
496  print*," APPLYING SOIL INCREMENTS FROM GSI OR JEDI"
497  ALLOCATE(stc_bck(lensfc, lsoil), smc_bck(lensfc, lsoil), slc_bck(lensfc,lsoil))
498  ALLOCATE(landinc_mask_fg(lensfc))
499  ENDIF
500  ! FOR NOW, CODE SO CAN DO BOTH, BUT MIGHT NEED TO THINK ABOUT THIS SOME MORE.
501  IF (do_sno_inc_jedi) THEN
502  ! ideally, would check here that sfcsub snow DA update is not also requested
503  ! but latter is controlled by fnsol, which is read in within that routine.
504  ! should be done at script level.
505  print*
506  print*," APPLYING SNOW INCREMENTS FROM JEDI"
507  ALLOCATE(snd_bck(lensfc), snd_inc(lensfc), swe_bck(lensfc))
508  ENDIF
509  ! set-up land mask info
510  ALLOCATE(landinc_mask(lensfc))
511  if (ivegsrc == 2) then ! sib
512  veg_type_landice=13
513  else
514  veg_type_landice=15
515  endif
516 ENDIF
517 
518 !--------------------------------------------------------------------------------
519 ! READ THE INPUT SURFACE DATA ON THE CUBED-SPHERE TILE.
520 !--------------------------------------------------------------------------------
521 
522  CALL read_data(lsoil,lensfc,do_nsst,do_sno_inc_jedi,do_soi_inc_jedi,.false.,is_noahmp=is_noahmp, &
523  tsffcs=tsffcs,smcfcs=smcfcs, &
524  swefcs=swefcs,stcfcs=stcfcs,tg3fcs=tg3fcs,zorfcs=zorfcs, &
525  cvfcs=cvfcs, cvbfcs=cvbfcs,cvtfcs=cvtfcs,albfcs=albfcs, &
526  vegfcs=vegfcs,slifcs=slifcs,cnpfcs=cnpfcs,f10m=f10m , &
527  vetfcs=vetfcs,sotfcs=sotfcs,alffcs=alffcs,ustar=ustar , &
528  fmm=fmm ,fhh=fhh ,sihfcs=sihfcs,sicfcs=sicfcs, &
529  sitfcs=sitfcs,tprcp=tprcp ,srflag=srflag,sndfcs=sndfcs, &
530  vmnfcs=vmnfcs,vmxfcs=vmxfcs,slcfcs=slcfcs,slpfcs=slpfcs, &
531  absfcs=absfcs,t2m=t2m ,q2m=q2m ,slmask=slmask, &
532  zsoil=zsoil, nsst=nsst)
533 
534  IF (frac_grid .AND. .NOT. is_noahmp) THEN
535  print *, 'FATAL ERROR: NOAH lsm update does not work with fractional grids.'
536  call mpi_abort(mpi_comm_world, 18, ierr)
537  ENDIF
538 
539  IF (is_noahmp .AND. do_sno_inc_jedi) THEN
540  print *, 'FATAL ERROR: Snow increment update does not work with NOAH_MP.'
541  call mpi_abort(mpi_comm_world, 29, ierr)
542  ENDIF
543 
544  IF (is_noahmp) THEN
545  lsm=lsm_noahmp
546  ELSE
547  lsm=lsm_noah
548  ENDIF
549 
550  IF (use_ufo) THEN
551  print*
552  print*,'USE UNFILTERED OROGRAPHY.'
553  ELSE
554  orog_uf = 0.0
555  ENDIF
556 
557  DO i=1,lensfc
558  aisfcs(i) = 0.
559  IF(nint(slifcs(i)).EQ.2) aisfcs(i) = 1.
560  ENDDO
561 
562  IF (do_nsst) THEN
563  sicfcs_fg=sicfcs
564  IF (.NOT. do_sfccycle ) THEN
565  print*
566  print*,"FIRST GUESS MASK ADJUSTED BY IFD RECORD"
567  slifcs_fg = slifcs
568  WHERE(nint(nsst%IFD) == 3) slifcs_fg = 2.0
569  ELSE
570  print*
571  print*,"SAVE FIRST GUESS MASK"
572  slifcs_fg = slifcs
573  ENDIF
574  ENDIF
575 
576  ! CALCULATE MASK FOR LAND INCREMENTS
577  IF (do_lndinc) &
578  CALL calculate_landinc_mask(swefcs, vetfcs, sotfcs, &
579  lensfc,veg_type_landice, landinc_mask)
580 
581 !--------------------------------------------------------------------------------
582 ! UPDATE SURFACE FIELDS.
583 !
584 ! FIRST, SET WATER AND LAND MASKS - SLMASKW/SLMASKL. FOR UNCOUPLED
585 ! (NON-FRACTIONAL) MODE, THESE ARE IDENTICAL TO THE CURRENT
586 ! MASK - '0' WATER; '1' LAND.
587 !--------------------------------------------------------------------------------
588 
589  IF (do_sfccycle) THEN
590  ALLOCATE(slmaskl(lensfc), slmaskw(lensfc))
591 
592  set_mask : IF (frac_grid) THEN
593 
594  DO i=1,lensfc
595  IF(landfrac(i) > 0.0_kind_io8) THEN
596  slmaskl(i) = ceiling(landfrac(i)-1.0e-6_kind_io8)
597  slmaskw(i) = floor(landfrac(i)+1.0e-6_kind_io8)
598  ELSE
599  IF(nint(slmask(i)) == 1) THEN ! If landfrac is zero, this should not happen.
600  ! So, stop processing.
601  print*, 'FATAL ERROR: LAND FRAC AND SLMASK MISMATCH.'
602  CALL mpi_abort(mpi_comm_world, 27, ierr)
603  ELSE
604  slmaskl(i) = 0.0_kind_io8
605  slmaskw(i) = 0.0_kind_io8
606  ENDIF
607  ENDIF
608 
609  ENDDO
610 
611  ELSE
612 
613 ! For running uncoupled (non-fractional grid).
614 
615  DO i=1,lensfc
616  IF(nint(slmask(i)) == 1) THEN
617  slmaskl(i) = 1.0_kind_io8
618  slmaskw(i) = 1.0_kind_io8
619  ELSE
620  slmaskl(i) = 0.0_kind_io8
621  slmaskw(i) = 0.0_kind_io8
622  ENDIF
623  ENDDO
624 
625  ENDIF set_mask
626 
627  DO i=1,lensfc
628  if(nint(slmask(i)) == 0) then
629  min_ice(i) = 0.15_kind_io8
630  else
631  min_ice(i) = 0.0_kind_io8
632  endif
633  ENDDO
634 
635  socfcs=0 ! Soil color. Not used yet.
636 
637  num_threads = num_parthds()
638  print*
639  print*,"CALL SFCCYCLE TO UPDATE SURFACE FIELDS."
640  CALL sfccycle(lugb,lensfc,lsoil,sig1t,deltsfc, &
641  iy,im,id,ih,fh,rla,rlo, &
642  slmaskl,slmaskw, orog, orog_uf, use_ufo, do_nsst, &
643  sihfcs,sicfcs,sitfcs,sndfcs,slcfcs, &
644  vmnfcs,vmxfcs,slpfcs,absfcs, &
645  tsffcs,swefcs,zorfcs,albfcs,tg3fcs, &
646  cnpfcs,smcfcs,stcfcs,slifcs,aisfcs, &
647  vegfcs,vetfcs,sotfcs,socfcs,alffcs, &
648  cvfcs,cvbfcs,cvtfcs,myrank,num_threads, nlunit, &
649  sz_nml, input_nml_file, &
650  min_ice, &
651  ialb,isot,ivegsrc,tile_num,i_index,j_index)
652 
653  DEALLOCATE(slmaskl, slmaskw)
654  ENDIF
655 
656 !--------------------------------------------------------------------------------
657 ! IF RUNNING WITH NSST, READ IN GSI FILE WITH THE UPDATED INCREMENTS (ON THE
658 ! GAUSSIAN GRID), INTERPOLATE INCREMENTS TO THE CUBED-SPHERE TILE, AND PERFORM
659 ! REQUIRED ADJUSTMENTS AND QC.
660 !--------------------------------------------------------------------------------
661 
662  IF (do_nsst) THEN
663  IF (nst_file == "NULL") THEN
664  print*
665  print*,"NO GSI FILE. ADJUST IFD FOR FORMER ICE POINTS."
666  DO i = 1, lensfc
667  IF (sicfcs_fg(i) > 0.0 .AND. sicfcs(i) == 0) THEN
668  nsst%IFD(i) = 3.0
669  ENDIF
670  ENDDO
671  nsst%TFINC = 0.0
672  ELSE
673  print*
674  print*,"ADJUST TREF FROM GSI INCREMENT"
675 !
676 ! Get tf climatology at the time
677 !
678  call get_tf_clm(rla,rlo,jdim,idim,iy,im,id,ih,tf_clm,tf_trd)
679  tf_clm_tile(:) = reshape(tf_clm, (/lensfc/) )
680  tf_trd_tile(:) = reshape(tf_trd, (/lensfc/) )
681 !
682 ! Get salinity climatology at the time
683 !
684  call get_sal_clm(rla,rlo,jdim,idim,iy,im,id,ih,sal_clm)
685  sal_clm_tile(:) = reshape(sal_clm, (/lensfc/) )
686 !
687 ! read tf analysis increment generated by GSI
688 !
689  CALL read_gsi_data(nst_file, 'NST')
690 !
691 ! update foundation & surface temperature for NSST
692 !
693  CALL adjust_nsst(rla,rlo,slifcs,slifcs_fg,tsffcs,sitfcs,sicfcs,sicfcs_fg,&
694  stcfcs,nsst,lensfc,lsoil,idim,jdim,zsea1,zsea2, &
695  tf_clm_tile,tf_trd_tile,sal_clm_tile,landfrac,frac_grid)
696  ENDIF
697  ENDIF
698 
699 !--------------------------------------------------------------------------------
700 ! READ IN AND APPLY LAND INCREMENTS FROM THE GSI/JEDI
701 !--------------------------------------------------------------------------------
702 
703  IF (do_lndinc) THEN
704 
705  ! SNOW INCREMENTS
706  ! do snow first, as temperature updates will use snow analaysis
707  IF (do_sno_inc_jedi) THEN
708  ! updates are made to snow depth only over land (and not-land ice).
709  ! SWE is then updated from the snow depth analysis, using the model
710  ! forecast density
711 
712  ! make sure incr. files exist
713  WRITE(rankch, '(I3.3)') (myrank+1)
714  jedi_sno_file = "snow_xainc." // rankch
715 
716  INQUIRE(file=trim(jedi_sno_file), exist=file_exists)
717  IF (.not. file_exists) then
718  print *, 'FATAL ERROR: snow increment (fv3 grid) update requested, &
719  but file does not exist : ', trim(jedi_sno_file)
720  call mpi_abort(mpi_comm_world, 10, ierr)
721  ENDIF
722 
723  !--------------------------------------------------------------------------------
724  ! read increments in
725  !--------------------------------------------------------------------------------
726 
727  ! Only coded for DA on native model grid (would always be the case for cycling DA)
728  CALL read_data(lsoil,lensfc,.false.,.true.,.false.,.true.,sndfcs=snd_inc)
729 
730  !--------------------------------------------------------------------------------
731  ! add increments to state vars
732  !--------------------------------------------------------------------------------
733 
734  ! store background states
735  snd_bck = sndfcs
736  swe_bck = swefcs
737 
738  CALL add_jedi_increment_snow(snd_inc,landinc_mask,lensfc,sndfcs)
739 
740  !--------------------------------------------------------------------------------
741  ! make any necessary adjustments to dependent variables
742  !--------------------------------------------------------------------------------
743 
744  CALL apply_land_da_adjustments_snd(lsm, lensfc, landinc_mask, swe_bck, snd_bck, &
745  sndfcs, swefcs)
746 
747  ENDIF ! jedi snow increments
748 
749  !re-calculate soilsnow mask if snow has been updated.
750  landinc_mask_fg = landinc_mask
751 
752  IF (do_sfccycle .OR. do_sno_inc_jedi) THEN
753  CALL calculate_landinc_mask(swefcs, vetfcs, sotfcs, lensfc, &
754  veg_type_landice, landinc_mask)
755  ENDIF
756 
757  ! store background states
758  stc_bck = stcfcs
759  smc_bck = smcfcs
760  slc_bck = slcfcs
761 
762  ! SOIL INCREMENTS
763  IF (do_soi_inc_gsi) THEN
764 
765  !--------------------------------------------------------------------------------
766  ! read increments in
767  !--------------------------------------------------------------------------------
768  ! make sure incr. files exist
769  WRITE(rankch, '(I3.3)') (myrank+1)
770  gsi_soi_file = "sfcincr_gsi." // rankch
771 
772  INQUIRE(file=trim(gsi_soi_file), exist=file_exists)
773  IF (.not. file_exists) then
774  print *, 'FATAL ERROR: gsi soil increment (gaussian grid) update requested, &
775  but file does not exist : ', trim(gsi_soi_file)
776  call mpi_abort(mpi_comm_world, 10, ierr)
777  ENDIF
778 
779  CALL read_gsi_data(gsi_soi_file, 'LND', lsoil=lsoil)
780 
781  !--------------------------------------------------------------------------------
782  ! interpolate increments to cubed sphere tiles
783  !--------------------------------------------------------------------------------
784 
785  CALL gaussian_to_fv3_interp(lsoil_incr,rla,rlo,&
786  stcinc,slcinc,landinc_mask,lensfc,lsoil,idim,jdim,lsm,myrank)
787 
788  !--------------------------------------------------------------------------------
789  ! save interpolated increments
790  !--------------------------------------------------------------------------------
791  CALL write_data(lensfc,idim,jdim,lsoil,do_nsst,.true.,nsst, &
792  stcinc=stcinc,slcinc=slcinc)
793 
794  ENDIF ! end reading and interpolating gsi soil increments
795 
796  IF (do_soi_inc_jedi) THEN
797 
798  !--------------------------------------------------------------------------------
799  ! read increments in
800  !--------------------------------------------------------------------------------
801  ! make sure incr. files exist
802  WRITE(rankch, '(I3.3)') (myrank+1)
803  jedi_soi_file = "soil_xainc." // rankch
804 
805  INQUIRE(file=trim(jedi_soi_file), exist=file_exists)
806  IF (.not. file_exists) then
807  print *, 'FATAL ERROR: soil increment (fv3 grid) update requested, but file &
808  does not exist: ', trim(jedi_soi_file)
809  call mpi_abort(mpi_comm_world, 10, ierr)
810  ENDIF
811 
812  CALL read_data(lsoil,lensfc,.false.,.false.,.true., &
813  .true.,is_noahmp=is_noahmp, &
814  stcinc=stcinc,slcinc=slcinc)
815 
816  ENDIF ! end reading jedi soil increments
817 
818  IF (do_soi_inc_gsi .or. do_soi_inc_jedi) THEN
819 
820  !--------------------------------------------------------------------------------
821  ! add increments to state vars
822  !--------------------------------------------------------------------------------
823 
824  ! below updates [STC/SMC/STC]FCS to hold the analysis
825  CALL add_increment_soil(lsoil_incr,stcinc,slcinc,stcfcs,smcfcs,slcfcs,stc_updated, &
826  slc_updated,landinc_mask,landinc_mask_fg,lensfc,lsoil,lsm,myrank)
827 
828  !--------------------------------------------------------------------------------
829  ! make any necessary adjustments to dependent variables
830  !--------------------------------------------------------------------------------
831 
832  CALL apply_land_da_adjustments_soil(lsoil_incr, lsm, isot, ivegsrc,lensfc, lsoil, &
833  sotfcs, landinc_mask_fg, stc_bck, stcfcs, smcfcs, slcfcs, stc_updated, &
834  slc_updated,zsoil)
835 
836 
837  ENDIF ! end applying soil increments and making adjustments
838 
839 !--------------------------------------------------------------------------------
840 ! clean up
841 !--------------------------------------------------------------------------------
842 
843  IF(ALLOCATED(landinc_mask_fg)) DEALLOCATE(landinc_mask_fg)
844  IF(ALLOCATED(landinc_mask)) DEALLOCATE(landinc_mask)
845  IF(ALLOCATED(stc_bck)) DEALLOCATE(stc_bck)
846  IF(ALLOCATED(smc_bck)) DEALLOCATE(smc_bck)
847  IF(ALLOCATED(slc_bck)) DEALLOCATE(slc_bck)
848  IF(ALLOCATED(snd_bck)) DEALLOCATE(snd_bck)
849  IF(ALLOCATED(swe_bck)) DEALLOCATE(swe_bck)
850  IF(ALLOCATED(snd_inc)) DEALLOCATE(snd_inc)
851 
852  ENDIF
853 !--------------------------------------------------------------------------------
854 ! WRITE OUT UPDATED SURFACE DATA ON THE CUBED-SPHERE TILE.
855 !--------------------------------------------------------------------------------
856 
857  IF (lsm==lsm_noahmp) THEN
858 
859  CALL write_data(lensfc,idim,jdim,lsoil,do_nsst,.false.,nsst,vegfcs=vegfcs, &
860  slcfcs=slcfcs,smcfcs=smcfcs,stcfcs=stcfcs,&
861  sicfcs=sicfcs,sihfcs=sihfcs)
862 
863  ELSEIF (lsm==lsm_noah) THEN
864 
865  CALL write_data(lensfc,idim,jdim,lsoil, &
866  do_nsst,.false.,nsst,slifcs=slifcs,tsffcs=tsffcs,vegfcs=vegfcs, &
867  swefcs=swefcs,tg3fcs=tg3fcs,zorfcs=zorfcs, &
868  albfcs=albfcs,alffcs=alffcs,cnpfcs=cnpfcs, &
869  f10m=f10m,t2m=t2m,q2m=q2m,vetfcs=vetfcs, &
870  sotfcs=sotfcs,ustar=ustar,fmm=fmm,fhh=fhh, &
871  sicfcs=sicfcs,sihfcs=sihfcs,sitfcs=sitfcs,tprcp=tprcp, &
872  srflag=srflag,swdfcs=sndfcs,vmnfcs=vmnfcs, &
873  vmxfcs=vmxfcs,slpfcs=slpfcs,absfcs=absfcs, &
874  slcfcs=slcfcs,smcfcs=smcfcs,stcfcs=stcfcs)
875 
876  ENDIF
877 
878  IF (do_nsst) THEN
879  DEALLOCATE(nsst%C_0)
880  DEALLOCATE(nsst%C_D)
881  DEALLOCATE(nsst%D_CONV)
882  DEALLOCATE(nsst%DT_COOL)
883  DEALLOCATE(nsst%IFD)
884  DEALLOCATE(nsst%QRAIN)
885  DEALLOCATE(nsst%TREF)
886  DEALLOCATE(nsst%TFINC)
887  DEALLOCATE(nsst%W_0)
888  DEALLOCATE(nsst%W_D)
889  DEALLOCATE(nsst%XS)
890  DEALLOCATE(nsst%XT)
891  DEALLOCATE(nsst%XTTS)
892  DEALLOCATE(nsst%XU)
893  DEALLOCATE(nsst%XV)
894  DEALLOCATE(nsst%XZ)
895  DEALLOCATE(nsst%XZTS)
896  DEALLOCATE(nsst%Z_C)
897  DEALLOCATE(nsst%ZM)
898  DEALLOCATE(slifcs_fg)
899  DEALLOCATE(sicfcs_fg)
900  ENDIF
901 
902  RETURN
903 
904  END SUBROUTINE sfcdrv
905 
937  SUBROUTINE adjust_nsst(RLA,RLO,SLMSK_TILE,SLMSK_FG_TILE,SKINT_TILE,&
938  SICET_TILE,sice_tile,sice_fg_tile,SOILT_TILE,NSST, &
939  LENSFC,LSOIL,IDIM,JDIM,ZSEA1,ZSEA2, &
940  tf_clm_tile,tf_trd_tile,sal_clm_tile,LANDFRAC, &
941  FRAC_GRID)
943  USE utils
944  USE gdswzd_mod
945  USE read_write_data, ONLY : idim_gaus, jdim_gaus, &
947  nsst_data
948 
949  USE mpi
950 
951  IMPLICIT NONE
952 
953  INTEGER, INTENT(IN) :: LENSFC, LSOIL, IDIM, JDIM
954 
955  LOGICAL, INTENT(IN) :: FRAC_GRID
956 
957  REAL, INTENT(IN) :: SLMSK_TILE(LENSFC), SLMSK_FG_TILE(LENSFC), LANDFRAC(LENSFC)
958  real, intent(in) :: tf_clm_tile(lensfc),tf_trd_tile(lensfc),sal_clm_tile(lensfc)
959  REAL, INTENT(IN) :: ZSEA1, ZSEA2,sice_tile(lensfc),sice_fg_tile(lensfc)
960  REAL, INTENT(IN) :: RLA(LENSFC), RLO(LENSFC)
961  REAL, INTENT(INOUT) :: SKINT_TILE(LENSFC)
962  REAL, INTENT(INOUT) :: SICET_TILE(LENSFC),SOILT_TILE(LENSFC,LSOIL)
963 
964  TYPE(NSST_DATA) :: NSST
965 
966  REAL, PARAMETER :: TMAX=313.0,tzero=273.16
967 
968  INTEGER :: IOPT, NRET, KGDS_GAUS(200)
969  INTEGER :: IGAUS, JGAUS, IJ, II, JJ, III, JJJ, KRAD
970  INTEGER :: ISTART, IEND, JSTART, JEND
971 !INTEGER :: MASK_TILE, MASK_FG_TILE
972  INTEGER,allocatable :: MASK_TILE(:),MASK_FG_TILE(:)
973  INTEGER :: ITILE, JTILE
974  INTEGER :: MAX_SEARCH, J, IERR
975  INTEGER :: IGAUSP1, JGAUSP1
976  integer :: nintp,nsearched,nice,nland
977  integer :: nfill,nfill_tice,nfill_clm
978  integer :: nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c
979 
980  INTEGER, ALLOCATABLE :: ID1(:,:), ID2(:,:), JDC(:,:)
981 
982  LOGICAL :: IS_ICE
983 
984  real :: tfreez
985  REAL :: TREF_SAVE,WSUM,tf_ice,tf_thaw
986  REAL :: FILL, DTZM, GAUS_RES_KM, DTREF
987  REAL, ALLOCATABLE :: XPTS(:), YPTS(:), LATS(:), LONS(:)
988  REAL, ALLOCATABLE :: DUM2D(:,:), LATS_RAD(:), LONS_RAD(:)
989  REAL, ALLOCATABLE :: AGRID(:,:,:), S2C(:,:,:)
990 
991  kgds_gaus = 0
992  kgds_gaus(1) = 4 ! OCT 6 - TYPE OF GRID (GAUSSIAN)
993  kgds_gaus(2) = idim_gaus ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
994  kgds_gaus(3) = jdim_gaus
995  kgds_gaus(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
996  kgds_gaus(5) = 0 ! OCT 14-16 - LON OF ORIGIN
997  kgds_gaus(6) = 128 ! OCT 17 - RESOLUTION FLAG
998  kgds_gaus(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
999  kgds_gaus(8) = nint(-360000./float(idim_gaus)) ! OCT 21-23 - LON OF EXTREME POINT
1000  kgds_gaus(9) = nint((360.0 / float(idim_gaus))*1000.0)
1001  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
1002  kgds_gaus(10) = jdim_gaus/2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
1003  kgds_gaus(12) = 255 ! OCT 29 - RESERVED
1004  kgds_gaus(20) = 255 ! OCT 5 - NOT USED, SET TO 255
1005 
1006  print*
1007  print*,'ADJUST NSST USING GSI INCREMENTS ON GAUSSIAN GRID'
1008 
1009 !----------------------------------------------------------------------
1010 ! CALL GDSWZD TO COMPUTE THE LAT/LON OF EACH GSI GAUSSIAN GRID POINT.
1011 !----------------------------------------------------------------------
1012 
1013  iopt = 0
1014  fill = -9999.
1015  ALLOCATE(xpts(idim_gaus*jdim_gaus))
1016  ALLOCATE(ypts(idim_gaus*jdim_gaus))
1017  ALLOCATE(lats(idim_gaus*jdim_gaus))
1018  ALLOCATE(lons(idim_gaus*jdim_gaus))
1019  xpts = fill
1020  ypts = fill
1021  lats = fill
1022  lons = fill
1023 
1024  CALL gdswzd(kgds_gaus,iopt,(idim_gaus*jdim_gaus),fill,xpts,ypts,lons,lats,nret)
1025 
1026  IF (nret /= (idim_gaus*jdim_gaus)) THEN
1027  print*,'FATAL ERROR: PROBLEM IN GDSWZD. STOP.'
1028  CALL mpi_abort(mpi_comm_world, 12, ierr)
1029  ENDIF
1030 
1031  DEALLOCATE (xpts, ypts)
1032 
1033  ALLOCATE(dum2d(idim_gaus,jdim_gaus))
1034  dum2d = reshape(lats, (/idim_gaus,jdim_gaus/) )
1035  DEALLOCATE(lats)
1036 
1037  ALLOCATE(lats_rad(jdim_gaus))
1038  DO j = 1, jdim_gaus
1039  lats_rad(j) = dum2d(1,jdim_gaus-j+1) * 3.1415926 / 180.0
1040  ENDDO
1041 
1042  dum2d = reshape(lons, (/idim_gaus,jdim_gaus/) )
1043  DEALLOCATE(lons)
1044  ALLOCATE(lons_rad(idim_gaus))
1045  lons_rad = dum2d(:,1) * 3.1415926 / 180.0
1046 
1047  DEALLOCATE(dum2d)
1048 
1049  ALLOCATE(agrid(idim,jdim,2))
1050  agrid(:,:,1) = reshape(rlo, (/idim,jdim/) )
1051  agrid(:,:,2) = reshape(rla, (/idim,jdim/) )
1052  agrid = agrid * 3.1415926 / 180.0
1053 
1054  ALLOCATE(id1(idim,jdim))
1055  ALLOCATE(id2(idim,jdim))
1056  ALLOCATE(jdc(idim,jdim))
1057  ALLOCATE(s2c(idim,jdim,4))
1058 
1059 !----------------------------------------------------------------------
1060 ! COMPUTE BILINEAR WEIGHTS FOR EACH MODEL POINT FROM THE NEAREST
1061 ! FOUR GSI/GAUSSIAN POINTS. DOES NOT ACCOUNT FOR MASK. THAT
1062 ! HAPPENS LATER.
1063 !----------------------------------------------------------------------
1064 
1065  CALL remap_coef( 1, idim, 1, jdim, idim_gaus, jdim_gaus, &
1066  lons_rad, lats_rad, id1, id2, jdc, s2c, agrid )
1067 
1068  DEALLOCATE(lons_rad, lats_rad, agrid)
1069 
1070 !----------------------------------------------------------------------
1071 ! THE MAXIMUM DISTANCE TO SEARCH IS 500 KM. HOW MANY GAUSSIAN
1072 ! GRID LENGTHS IS THAT?
1073 !----------------------------------------------------------------------
1074 
1075  gaus_res_km = 360.0 / idim_gaus * 111.0
1076  max_search = ceiling(500.0/gaus_res_km)
1077 
1078  print*
1079  print*,'MAXIMUM SEARCH IS ',max_search, ' GAUSSIAN POINTS.'
1080  print*
1081 
1082 !
1083 ! Initialize variables for counts statitics to be zeros
1084 !
1085  nintp = 0
1086  nset_thaw = 0
1087  nset_thaw_s = 0
1088  nset_thaw_i = 0
1089  nset_thaw_c = 0
1090  nsearched = 0
1091  nfill = 0
1092  nfill_tice = 0
1093  nfill_clm = 0
1094  nice = 0
1095  nland = 0
1096 !----------------------------------------------------------------------
1097 ! TREF INCREMENT WILL BE OUTPUT. INITIALIZE TO ZERO.
1098 !----------------------------------------------------------------------
1099 
1100  nsst%TFINC = 0.0
1101 
1102  allocate(mask_tile(lensfc))
1103  allocate(mask_fg_tile(lensfc))
1104 
1105  IF(.NOT. frac_grid) THEN
1106  mask_tile = nint(slmsk_tile)
1107  mask_fg_tile = nint(slmsk_fg_tile)
1108  ELSE
1109  mask_tile=0
1110  WHERE(sice_tile > 0.0) mask_tile=2
1111  WHERE(landfrac == 1.0) mask_tile=1
1112  mask_fg_tile=0
1113  WHERE(sice_fg_tile > 0.0) mask_fg_tile=2
1114  WHERE(landfrac == 1.0) mask_fg_tile=1
1115  ENDIF
1116 
1117  ij_loop : DO ij = 1, lensfc
1118 
1119 !
1120 ! when sea ice exists, get salinity dependent water temperature
1121 !
1122  tf_ice = tfreez(sal_clm_tile(ij)) + tzero
1123 !----------------------------------------------------------------------
1124 ! SKIP LAND POINTS. NSST NOT APPLIED AT LAND.
1125 !----------------------------------------------------------------------
1126 
1127  IF (mask_tile(ij) == 1) THEN
1128  nland = nland + 1
1129  cycle ij_loop
1130  ENDIF
1131 
1132 !
1133 ! these are ice points. set tref to tf_ice and update tmpsfc.
1134 !
1135  if (mask_tile(ij) == 2) then
1136  nsst%tref(ij)=tf_ice ! water part tmp set
1137  skint_tile(ij)=(1.0-sice_tile(ij))*nsst%tref(ij)+sice_tile(ij)*sicet_tile(ij)
1138  nice = nice + 1
1139  cycle ij_loop
1140  endif
1141 
1142 !
1143 ! Get i,j index on array of (idim,jdim) from known ij
1144 !
1145  jtile = (ij-1) / idim + 1
1146  itile = mod(ij,idim)
1147  IF (itile==0) itile = idim
1148 
1149 !----------------------------------------------------------------------
1150 ! IF THE MODEL POINT WAS ICE COVERED, BUT IS NOW OPEN WATER, SET
1151 ! TREF TO searched adjascent open water onea, if failed the search, set to
1152 ! weighted average of tf_ice and tf_clm. For NSST vars, set xz TO '30' AND ALL OTHER FIELDS TO ZERO.
1153 !----------------------------------------------------------------------
1154 
1155  IF (mask_fg_tile(ij) == 2 .AND. mask_tile(ij) == 0) THEN
1156 !
1157 ! set background for the thaw (just melted water) situation
1158 !
1159  call tf_thaw_set(nsst%tref,mask_fg_tile,itile,jtile,tf_ice,tf_clm_tile(ij),tf_thaw,idim,jdim, &
1160  nset_thaw_s,nset_thaw_i,nset_thaw_c)
1161  call nsst_water_reset(nsst,ij,tf_thaw)
1162  nset_thaw = nset_thaw + 1
1163  ENDIF
1164 
1165 !----------------------------------------------------------------------
1166 ! THESE ARE POINTS THAT ARE OPEN WATER AND WERE OPEN WATER PRIOR
1167 ! TO ANY ICE UPDATE BY SFCCYCLE. UPDATE TREF AND SKIN TEMP.
1168 ! AT OPEN WATER POINTS, THE SEA ICE TEMPERATURE (SICET_TILE) AND
1169 ! SOIL COLUMN TEMPERATURE (SOILT_TILE) ARE SET TO THE SKIN TEMP.
1170 ! IT IS SIMPLY A FILLER VALUE. THESE FIELDS ARE NOT USED AT
1171 ! OPEN WATER POINTS.
1172 !----------------------------------------------------------------------
1173 !----------------------------------------------------------------------
1174 ! SEE IF ANY OF THE NEAREST GSI POINTS MASK AREA OPEN WATER.
1175 ! IF SO, APPLY NSST INCREMENT USING BILINEAR INTERPOLATION.
1176 !----------------------------------------------------------------------
1177 
1178  igaus = id1(itile,jtile)
1179  jgaus = jdc(itile,jtile)
1180  igausp1 = id2(itile,jtile)
1181  jgausp1 = jdc(itile,jtile)+1
1182 
1183  IF (slmsk_gaus(igaus,jgaus) == 0 .OR. &
1184  slmsk_gaus(igausp1,jgaus) == 0 .OR. &
1185  slmsk_gaus(igausp1,jgausp1) == 0 .OR. &
1186  slmsk_gaus(igaus,jgausp1) == 0) THEN
1187 
1188  dtref = 0.0
1189  wsum = 0.0
1190 
1191  IF (slmsk_gaus(igaus,jgaus) == 0) THEN
1192  dtref = dtref + (s2c(itile,jtile,1) * dtref_gaus(igaus,jgaus))
1193  wsum = wsum + s2c(itile,jtile,1)
1194  ENDIF
1195 
1196  IF (slmsk_gaus(igausp1,jgaus) == 0) THEN
1197  dtref = dtref + (s2c(itile,jtile,2) * dtref_gaus(igausp1,jgaus))
1198  wsum = wsum + s2c(itile,jtile,2)
1199  ENDIF
1200 
1201  IF (slmsk_gaus(igausp1,jgausp1) == 0) THEN
1202  dtref = dtref + (s2c(itile,jtile,3) * dtref_gaus(igausp1,jgausp1))
1203  wsum = wsum + s2c(itile,jtile,3)
1204  ENDIF
1205 
1206  IF (slmsk_gaus(igaus,jgausp1) == 0) THEN
1207  dtref = dtref + (s2c(itile,jtile,4) * dtref_gaus(igaus,jgausp1))
1208  wsum = wsum + s2c(itile,jtile,4)
1209  ENDIF
1210 
1211  nintp = nintp + 1
1212  dtref = dtref / wsum
1213 
1214  tref_save = nsst%TREF(ij)
1215  nsst%TREF(ij) = nsst%TREF(ij) + dtref
1216  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1217  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1218  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1219 
1220  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1221  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1222 
1223  skint_tile(ij) = nsst%TREF(ij) + dtzm
1224  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1225  skint_tile(ij) = min(skint_tile(ij), tmax)
1226 
1227  sicet_tile(ij) = skint_tile(ij)
1228 ! Under fractional grids, soilt is used at points with at
1229 ! least some land.
1230  IF(.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1231 
1232 !----------------------------------------------------------------------
1233 ! NO NEARBY GSI/GAUSSIAN OPEN WATER POINTS. PERFORM A SPIRAL SEARCH TO
1234 ! FIND NEAREST NON-LAND POINT ON GSI/GAUSSIAN GRID.
1235 !----------------------------------------------------------------------
1236 
1237  ELSE
1238 
1239  is_ice = .false.
1240 
1241  DO krad = 1, max_search
1242 
1243  istart = igaus - krad
1244  iend = igaus + krad
1245  jstart = jgaus - krad
1246  jend = jgaus + krad
1247 
1248  DO jj = jstart, jend
1249  DO ii = istart, iend
1250 
1251  IF((jj == jstart) .OR. (jj == jend) .OR. &
1252  (ii == istart) .OR. (ii == iend)) THEN
1253 
1254  IF ((jj >= 1) .AND. (jj <= jdim_gaus)) THEN
1255 
1256  jjj = jj
1257  IF (ii <= 0) THEN
1258  iii = idim_gaus + ii
1259  ELSE IF (ii >= (idim_gaus+1)) THEN
1260  iii = ii - idim_gaus
1261  ELSE
1262  iii = ii
1263  END IF
1264 
1265 !----------------------------------------------------------------------
1266 ! SEE IF NEARBY POINTS ARE SEA ICE. IF THEY ARE, AND THE SEARCH FOR
1267 ! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
1268 ! FREEZING BELOW.
1269 !----------------------------------------------------------------------
1270 
1271  IF (krad <= 2 .AND. slmsk_gaus(iii,jjj) == 2) is_ice = .true.
1272 
1273  IF (slmsk_gaus(iii,jjj) == 0) THEN
1274 
1275 ! PRINT*,'MISMATCH AT TILE POINT ',ITILE,JTILE
1276 ! PRINT*,'UPDATE TREF USING GSI INCREMENT AT ',III,JJJ,DTREF_GAUS(III,JJJ)
1277  nsearched = nsearched + 1
1278 
1279  tref_save = nsst%TREF(ij)
1280  nsst%TREF(ij ) = nsst%TREF(ij) + dtref_gaus(iii,jjj)
1281  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1282  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1283  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1284 
1285  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1286  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1287 
1288  skint_tile(ij) = nsst%TREF(ij) + dtzm
1289  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1290  skint_tile(ij) = min(skint_tile(ij), tmax)
1291 
1292  sicet_tile(ij) = skint_tile(ij)
1293  IF(.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1294  cycle ij_loop
1295 
1296  ENDIF ! GSI/Gaussian mask is open water
1297 
1298  ENDIF
1299 
1300  ENDIF
1301 
1302  ENDDO
1303  ENDDO
1304 
1305  ENDDO ! KRAD LOOP
1306 
1307 !----------------------------------------------------------------------
1308 ! THE SEARCH FAILED. IF THERE IS NEARBY ICE, SET TREF TO FREEZING.
1309 ! ELSE UPDATE TREF BASED ON THE ANNUAL SST CYCLE.
1310 !----------------------------------------------------------------------
1311 
1312 ! PRINT*,'WARNING !!!!!! SEARCH FAILED AT TILE POINT ',ITILE,JTILE
1313 
1314  nfill = nfill + 1
1315  IF (is_ice) THEN
1316  nsst%TREF(ij) = tf_ice
1317 ! PRINT*,"NEARBY ICE. SET TREF TO FREEZING"
1318  nfill_tice = nfill_tice + 1
1319  ELSE
1320  tref_save = nsst%TREF(ij)
1321  nsst%TREF(ij) = nsst%TREF(ij) + tf_trd_tile(ij)
1322  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1323  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1324  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1325 ! PRINT*,'UPDATE TREF FROM SST CLIMO ',DTREF
1326  nfill_clm = nfill_clm + 1
1327  ENDIF
1328 
1329  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1330  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1331 
1332  skint_tile(ij) = nsst%TREF(ij) + dtzm
1333  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1334  skint_tile(ij) = min(skint_tile(ij), tmax)
1335 
1336  sicet_tile(ij) = skint_tile(ij)
1337  IF (.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1338 
1339  ENDIF ! NEARBY GAUSSIAN POINTS ARE OPEN WATER?
1340 
1341  ENDDO ij_loop
1342 
1343  write(*,'(a)') 'statistics of grids number processed for tile : '
1344  write(*,'(a,I8)') ' nintp = ',nintp
1345  write(*,'(a,4I8)') 'nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c =',nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c
1346  write(*,'(a,I8)') ' nsearched = ',nsearched
1347  write(*,'(a,3I6)') ' nfill,nfill_tice,nfill_clm = ',nfill,nfill_tice,nfill_clm
1348  write(*,'(a,I8)') ' nice = ',nice
1349  write(*,'(a,I8)') ' nland = ',nland
1350 
1351  DEALLOCATE(id1, id2, jdc, s2c, mask_tile, mask_fg_tile)
1352 
1353  END SUBROUTINE adjust_nsst
1354 
1365  SUBROUTINE climo_trend(LATITUDE, MON, DAY, DELTSFC, DTREF)
1366  IMPLICIT NONE
1367 
1368  INTEGER, INTENT(IN) :: MON, DAY
1369 
1370  REAL, INTENT(IN) :: LATITUDE, DELTSFC
1371  REAL, INTENT(OUT) :: DTREF
1372 
1373  INTEGER :: NUM_DAYS(12), MON2, MON1
1374 
1375  REAL, TARGET :: SST_80_90(12)
1376  REAL, TARGET :: SST_70_80(12)
1377  REAL, TARGET :: SST_60_70(12)
1378  REAL, TARGET :: SST_50_60(12)
1379  REAL, TARGET :: SST_40_50(12)
1380  REAL, TARGET :: SST_30_40(12)
1381  REAL, TARGET :: SST_20_30(12)
1382  REAL, TARGET :: SST_10_20(12)
1383  REAL, TARGET :: SST_00_10(12)
1384  REAL, TARGET :: SST_M10_00(12)
1385  REAL, TARGET :: SST_M20_M10(12)
1386  REAL, TARGET :: SST_M30_M20(12)
1387  REAL, TARGET :: SST_M40_M30(12)
1388  REAL, TARGET :: SST_M50_M40(12)
1389  REAL, TARGET :: SST_M60_M50(12)
1390  REAL, TARGET :: SST_M70_M60(12)
1391  REAL, TARGET :: SST_M80_M70(12)
1392  REAL, TARGET :: SST_M90_M80(12)
1393 
1394  REAL, POINTER :: SST(:)
1395 
1396  DATA num_days /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
1397 
1398  DATA sst_80_90 /271.466, 271.458, 271.448, 271.445, 271.519, 271.636, &
1399  272.023, 272.066, 272.001, 271.698, 271.510, 271.472/
1400 
1401  DATA sst_70_80 /272.149, 272.103, 272.095, 272.126, 272.360, 272.988, &
1402  274.061, 274.868, 274.415, 273.201, 272.468, 272.268/
1403 
1404  DATA sst_60_70 /274.240, 274.019, 273.988, 274.185, 275.104, 276.875, &
1405  279.005, 280.172, 279.396, 277.586, 275.818, 274.803/
1406 
1407  DATA sst_50_60 /277.277, 276.935, 277.021, 277.531, 279.100, 281.357, &
1408  283.735, 285.171, 284.399, 282.328, 279.918, 278.199/
1409 
1410  DATA sst_40_50 /281.321, 280.721, 280.850, 281.820, 283.958, 286.588, &
1411  289.195, 290.873, 290.014, 287.652, 284.898, 282.735/
1412 
1413  DATA sst_30_40 /289.189, 288.519, 288.687, 289.648, 291.547, 293.904, &
1414  296.110, 297.319, 296.816, 295.225, 292.908, 290.743/
1415 
1416  DATA sst_20_30 /294.807, 294.348, 294.710, 295.714, 297.224, 298.703, &
1417  299.682, 300.127, 300.099, 299.455, 297.953, 296.177/
1418 
1419  DATA sst_10_20 /298.878, 298.720, 299.033, 299.707, 300.431, 300.709, &
1420  300.814, 300.976, 301.174, 301.145, 300.587, 299.694/
1421 
1422  DATA sst_00_10 /300.415, 300.548, 300.939, 301.365, 301.505, 301.141, &
1423  300.779, 300.660, 300.818, 300.994, 300.941, 300.675/
1424 
1425  DATA sst_m10_00 /300.226, 300.558, 300.914, 301.047, 300.645, 299.870, &
1426  299.114, 298.751, 298.875, 299.294, 299.721, 299.989/
1427 
1428  DATA sst_m20_m10 /299.547, 299.985, 300.056, 299.676, 298.841, 297.788, &
1429  296.893, 296.491, 296.687, 297.355, 298.220, 298.964/
1430 
1431  DATA sst_m30_m20 /297.524, 298.073, 297.897, 297.088, 295.846, 294.520, &
1432  293.525, 293.087, 293.217, 293.951, 295.047, 296.363/
1433 
1434  DATA sst_m40_m30 /293.054, 293.765, 293.468, 292.447, 291.128, 289.781, &
1435  288.773, 288.239, 288.203, 288.794, 289.947, 291.553/
1436 
1437  DATA sst_m50_m40 /285.052, 285.599, 285.426, 284.681, 283.761, 282.826, &
1438  282.138, 281.730, 281.659, 281.965, 282.768, 283.961/
1439 
1440  DATA sst_m60_m50 /277.818, 278.174, 277.991, 277.455, 276.824, 276.229, &
1441  275.817, 275.585, 275.560, 275.687, 276.142, 276.968/
1442 
1443  DATA sst_m70_m60 /273.436, 273.793, 273.451, 272.813, 272.349, 272.048, &
1444  271.901, 271.838, 271.845, 271.889, 272.080, 272.607/
1445 
1446  DATA sst_m80_m70 /271.579, 271.578, 271.471, 271.407, 271.392, 271.391, &
1447  271.390, 271.391, 271.394, 271.401, 271.422, 271.486/
1448 
1449  DATA sst_m90_m80 /271.350, 271.350, 271.350, 271.350, 271.350, 271.350, &
1450  271.350, 271.350, 271.350, 271.350, 271.350, 271.350/
1451 
1452  NULLIFY(sst)
1453  IF (latitude > 80.0) THEN
1454  sst => sst_80_90
1455  ELSEIF (latitude > 70.0) THEN
1456  sst => sst_70_80
1457  ELSEIF (latitude > 60.0) THEN
1458  sst => sst_60_70
1459  ELSEIF (latitude > 50.0) THEN
1460  sst => sst_50_60
1461  ELSEIF (latitude > 40.0) THEN
1462  sst => sst_40_50
1463  ELSEIF (latitude > 30.0) THEN
1464  sst => sst_30_40
1465  ELSEIF (latitude > 20.0) THEN
1466  sst => sst_20_30
1467  ELSEIF (latitude > 10.0) THEN
1468  sst => sst_10_20
1469  ELSEIF (latitude > 0.0) THEN
1470  sst => sst_00_10
1471  ELSEIF (latitude > -10.0) THEN
1472  sst => sst_m10_00
1473  ELSEIF (latitude > -20.0) THEN
1474  sst => sst_m20_m10
1475  ELSEIF (latitude > -30.0) THEN
1476  sst => sst_m30_m20
1477  ELSEIF (latitude > -40.0) THEN
1478  sst => sst_m40_m30
1479  ELSEIF (latitude > -50.0) THEN
1480  sst => sst_m50_m40
1481  ELSEIF (latitude > -60.0) THEN
1482  sst => sst_m60_m50
1483  ELSEIF (latitude > -70.0) THEN
1484  sst => sst_m70_m60
1485  ELSEIF (latitude > -80.0) THEN
1486  sst => sst_m80_m70
1487  ELSE
1488  sst => sst_m90_m80
1489  END IF
1490 
1491  IF (day >= 15) THEN
1492  mon2 = mon+1
1493  IF(mon2 == 13) mon2 = 1
1494  mon1 = mon
1495  dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
1496  ELSE
1497  mon1 = mon - 1
1498  IF (mon1 == 0) mon1=12
1499  mon2 = mon
1500  dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
1501  ENDIF
1502 
1503  dtref = dtref * (deltsfc / 24.0)
1504 
1505  END SUBROUTINE climo_trend
1506 
1518  SUBROUTINE dtzm_point(XT,XZ,DT_COOL,ZC,Z1,Z2,DTZM)
1519  implicit none
1520 
1521  real, intent(in) :: xt,xz,dt_cool,zc,z1,z2
1522  real, intent(out) :: dtzm
1523 
1524  real, parameter :: zero = 0.0
1525  real, parameter :: one = 1.0
1526  real, parameter :: half = 0.5
1527  real :: dt_warm,dtw,dtc
1528 !
1529 ! get the mean warming in the range of z=z1 to z=z2
1530 !
1531  dtw = zero
1532  if ( xt > zero ) then
1533  dt_warm = (xt+xt)/xz ! Tw(0)
1534  if ( z1 < z2) then
1535  if ( z2 < xz ) then
1536  dtw = dt_warm*(one-(z1+z2)/(xz+xz))
1537  elseif ( z1 < xz .and. z2 >= xz ) then
1538  dtw = half*(one-z1/xz)*dt_warm*(xz-z1)/(z2-z1)
1539  endif
1540  elseif ( z1 == z2 ) then
1541  if ( z1 < xz ) then
1542  dtw = dt_warm*(one-z1/xz)
1543  endif
1544  endif
1545  endif
1546 !
1547 ! get the mean cooling in the range of z=z1 to z=z2
1548 !
1549  dtc = zero
1550  if ( zc > zero ) then
1551  if ( z1 < z2) then
1552  if ( z2 < zc ) then
1553  dtc = dt_cool*(one-(z1+z2)/(zc+zc))
1554  elseif ( z1 < zc .and. z2 >= zc ) then
1555  dtc = half*(one-z1/zc)*dt_cool*(zc-z1)/(z2-z1)
1556  endif
1557  elseif ( z1 == z2 ) then
1558  if ( z1 < zc ) then
1559  dtc = dt_cool*(one-z1/zc)
1560  endif
1561  endif
1562  endif
1563 
1564 !
1565 ! get the mean T departure from Tf in the range of z=z1 to z=z2
1566 !
1567  dtzm = dtw - dtc
1568 
1569  END SUBROUTINE dtzm_point
1570 
1590  subroutine tf_thaw_set(tf_ij,mask_ij,itile,jtile,tice,tclm,tf_thaw,nx,ny, &
1591  nset_thaw_s,nset_thaw_i,nset_thaw_c)
1593  real, dimension(nx*ny), intent(in) :: tf_ij
1594  integer, dimension(nx*ny), intent(in) :: mask_ij
1595  real, intent(in) :: tice,tclm
1596  integer, intent(in) :: itile,jtile,nx,ny
1597  real, intent(out) :: tf_thaw
1598  integer, intent(inout) :: nset_thaw_s,nset_thaw_i,nset_thaw_c
1599 ! Local
1600  real, parameter :: bmiss = -999.0
1601  real, dimension(nx,ny) :: tf
1602  integer, dimension(nx,ny) :: mask
1603  integer :: krad,max_search,istart,iend,jstart,jend
1604  integer :: ii,jj,iii,jjj
1605  logical :: is_ice
1606 
1607  max_search = 2
1608 
1609  mask(:,:) = reshape(mask_ij,(/nx,ny/) )
1610  tf(:,:) = reshape(tf_ij,(/nx,ny/) )
1611 
1612  tf_thaw = bmiss
1613 
1614  do krad = 1, max_search
1615 
1616  istart = itile - krad
1617  iend = itile + krad
1618  jstart = jtile - krad
1619  jend = jtile + krad
1620 
1621  do jj = jstart, jend
1622  do ii = istart, iend
1623 
1624 
1625  if ((jj == jstart) .or. (jj == jend) .or. &
1626  (ii == istart) .or. (ii == iend)) then
1627 
1628  if ((jj >= 1) .and. (jj <= ny)) then
1629  jjj = jj
1630  if (ii <= 0) then
1631  iii = nx + ii
1632  else if (ii >= (nx+1)) then
1633  iii = ii - nx
1634  else
1635  iii = ii
1636  endif
1637 
1638 !----------------------------------------------------------------------
1639 ! SEE IF NEARBY POINTS ARE SEA ICE. IF THEY ARE, AND THE SEARCH FOR
1640 ! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
1641 ! FREEZING BELOW.
1642 !----------------------------------------------------------------------
1643  if (krad <= 2 .and. mask(iii,jjj) == 2) is_ice = .true.
1644 
1645  if (mask(iii,jjj) == 0) then
1646  tf_thaw = tf(iii,jjj)
1647  nset_thaw_s = nset_thaw_s + 1
1648  write(*,'(a,I4,2F9.3)') 'nset_thaw_s,tf(iii,jjj),tclm : ',nset_thaw_s,tf(iii,jjj),tclm
1649  go to 100
1650  endif ! tile mask is open water
1651 
1652  endif
1653  endif
1654  enddo
1655  enddo
1656  enddo ! krad loop
1657 
1658  100 continue
1659 
1660  if ( tf_thaw == bmiss ) then
1661  if (is_ice) then
1662  tf_thaw = tice
1663  nset_thaw_i = nset_thaw_i + 1
1664  write(*,'(a,I4,F9.3)') 'nset_thaw_i,tf_ice : ',nset_thaw_i,tice
1665  else
1666  tf_thaw = 0.8*tice+0.2*tclm
1667  nset_thaw_c = nset_thaw_c + 1
1668  write(*,'(a,I4,2F9.3)') 'nset_thaw_c,tf_ice,tclm : ',nset_thaw_c,tice,tclm
1669  endif
1670  endif
1671 
1672  end subroutine tf_thaw_set
1673 
1681  subroutine nsst_water_reset(nsst,ij,tf_thaw)
1683  implicit none
1684 
1685  integer, intent(in) :: ij
1686 
1687  real, intent(in) :: tf_thaw
1688 
1689  type(nsst_data), intent(inout) :: nsst
1690 
1691  nsst%c_0(ij) = 0.0
1692  nsst%c_d(ij) = 0.0
1693  nsst%d_conv(ij) = 0.0
1694  nsst%dt_cool(ij) = 0.0
1695  nsst%ifd(ij) = 0.0
1696  nsst%qrain(ij) = 0.0
1697  nsst%tref(ij) = tf_thaw
1698  nsst%w_0(ij) = 0.0
1699  nsst%w_d(ij) = 0.0
1700  nsst%xs(ij) = 0.0
1701  nsst%xt(ij) = 0.0
1702  nsst%xtts(ij) = 0.0
1703  nsst%xu(ij) = 0.0
1704  nsst%xv(ij) = 0.0
1705  nsst%xz(ij) = 30.0
1706  nsst%xzts(ij) = 0.0
1707  nsst%z_c(ij) = 0.0
1708  nsst%zm(ij) = 0.0
1709 
1710  end subroutine nsst_water_reset
1711 
1727 subroutine get_tf_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,tf_clm,tf_trd)
1729 
1730  implicit none
1731 
1732  real, dimension(nx*ny), intent(in) :: xlats_ij
1733  real, dimension(nx*ny), intent(in) :: xlons_ij
1734  real, dimension(nx,ny), intent(out) :: tf_clm
1735  real, dimension(nx,ny), intent(out) :: tf_trd
1736  integer, intent(in) :: iy,im,id,ih,nx,ny
1737 ! local declare
1738  real, allocatable, dimension(:,:) :: tf_clm0 ! sst climatology at the valid time (nxc,nyc)
1739  real, allocatable, dimension(:,:) :: tf_trd0 ! 6-hourly sst climatology tendency valid at atime (nxc,nyc)
1740  real, allocatable, dimension(:) :: cxlats ! latitudes of sst climatology
1741  real, allocatable, dimension(:) :: cxlons ! longitudes of sst climatology
1742 
1743  real, dimension(nx*ny) :: tf_clm_ij ! sst climatology at target grids (nx*ny)
1744  real, dimension(nx*ny) :: tf_trd_ij ! 6-hourly sst climatology tendency
1745  real :: wei1,wei2
1746  integer :: nxc,nyc,mon1,mon2
1747  character (len=6), parameter :: fin_tf_clm='sstclm' ! sst climatology file name
1748 !
1749 ! get which two months used and their weights from atime
1750 !
1751  call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1752 !
1753 ! get the dimensions of the sst climatology & allocate the related arrays
1754 !
1755  call get_tf_clm_dim(fin_tf_clm,nyc,nxc)
1756  allocate( tf_clm0(nxc,nyc),tf_trd0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
1757 !
1758 ! get tf_clm at the analysis time from monthly climatology & cxlats, cxlons
1759 !
1760  call get_tf_clm_ta(tf_clm0,tf_trd0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
1761 !
1762 ! get tf_clm (nx by ny lat/lon) valid at atime
1763 !
1764  if ( nx == nxc .and. ny == nyc ) then
1765  tf_clm(:,:) = tf_clm0(:,:)
1766  tf_trd(:,:) = tf_trd0(:,:)
1767 ! write(*,'(a,2f9.3)') 'same dimensions, tf_clm, min : ',minval(tf_clm),maxval(tf_clm)
1768  else
1769 ! write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
1770  call intp_tile(tf_clm0, cxlats, cxlons, nyc, nxc, &
1771  tf_clm_ij,xlats_ij,xlons_ij,ny, nx)
1772  call intp_tile(tf_trd0, cxlats, cxlons, nyc, nxc, &
1773  tf_trd_ij,xlats_ij,xlons_ij,ny, nx)
1774 ! write(*,'(a,2f9.3)') 'tf_clm0, min, max : ',minval(tf_clm0),maxval(tf_clm0)
1775 
1776  tf_clm(:,:) = reshape(tf_clm_ij, (/nx,ny/) )
1777  tf_trd(:,:) = reshape(tf_trd_ij, (/nx,ny/) )
1778  endif
1779 
1780 end subroutine get_tf_clm
1781 
1796 subroutine get_tf_clm_ta(tf_clm_ta,tf_clm_trend,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
1798  implicit none
1799 
1800 ! input
1801  integer, intent(in) :: nlat,nlon,mon1,mon2
1802  real, intent(in) :: wei1,wei2
1803 ! output
1804  real, dimension(nlon,nlat), intent(out) :: tf_clm_ta,tf_clm_trend
1805  real, dimension(nlat), intent(out) :: xlats
1806  real, dimension(nlon), intent(out) :: xlons
1807 
1808 !input/output data file names
1809  character (len=6), parameter :: fin_tf_clm='sstclm'
1810 
1811 ! local declare
1812  real, dimension(nlon,nlat) :: tf_clm1,tf_clm2
1813 
1814 !
1815 ! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
1816 !
1817  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm1,xlats,xlons,nlat,nlon,mon1)
1818  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm2,xlats,xlons,nlat,nlon,mon2)
1819 !
1820 ! tf_clim at the analysis time
1821 !
1822  tf_clm_ta(:,:) = wei1*tf_clm1(:,:)+wei2*tf_clm2(:,:)
1823 !
1824 ! tf_clim trend at the analysis time (month to month)
1825 !
1826  tf_clm_trend(:,:) = (tf_clm2(:,:)-tf_clm1(:,:))/120.0
1827 
1828  write(*,'(a,2f9.3)') 'tf_clm_ta, min, max : ',minval(tf_clm_ta),maxval(tf_clm_ta)
1829  write(*,'(a,2f9.3)') 'tf_clm_trend, min, max : ',minval(tf_clm_trend),maxval(tf_clm_trend)
1830  end subroutine get_tf_clm_ta
1831 
1844 subroutine get_sal_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,sal_clm)
1846  implicit none
1847 
1848  real, dimension(nx*ny), intent(in) :: xlats_ij !
1849  real, dimension(nx*ny), intent(in) :: xlons_ij !
1850  real, dimension(nx,ny), intent(out) :: sal_clm !
1851  integer, intent(in) :: iy,im,id,ih,nx,ny
1852 ! local declare
1853  real, allocatable, dimension(:,:) :: sal_clm0 ! salinity climatology at the valid time
1854  real, allocatable, dimension(:) :: cxlats ! latitudes of sst climatology
1855  real, allocatable, dimension(:) :: cxlons ! longitudes of sst climatology
1856 
1857  real, dimension(nx*ny) :: sal_clm_ij ! salinity climatology at target grids (nx*ny)
1858  real :: wei1,wei2
1859  integer :: nxc,nyc,mon1,mon2
1860  character (len=6), parameter :: fin_sal_clm='salclm' ! salinity climatology file name
1861 !
1862 ! get which two months used and their weights from atime
1863 !
1864  call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1865 !
1866 ! get the dimensions of the sst climatology & allocate the related arrays
1867 !
1868  call get_dim_nc(fin_sal_clm,nyc,nxc)
1869  allocate( sal_clm0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
1870 !
1871 ! get sal_clm at the analysis time from monthly climatology & cxlats, cxlons
1872 !
1873  call get_sal_clm_ta(sal_clm0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
1874 !
1875 ! get sal_clm (nx by ny lat/lon) valid at atime
1876 !
1877  if ( nx == nxc .and. ny == nyc ) then
1878  sal_clm(:,:) = sal_clm0(:,:)
1879 ! write(*,'(a,2f9.3)') 'same dimensions, sal_clm, min : ',minval(sal_clm),maxval(sal_clm)
1880  else
1881 ! write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
1882  call intp_tile(sal_clm0, cxlats, cxlons, nyc, nxc, &
1883  sal_clm_ij,xlats_ij,xlons_ij,ny, nx)
1884 ! write(*,'(a,2f9.3)') 'sal_clm0, min, max : ',minval(sal_clm0),maxval(sal_clm0)
1885 ! write(*,'(a,2f9.3)') 'done with intp_tile for sal_clm, min, max : ',minval(sal_clm_ij),maxval(sal_clm_ij)
1886 
1887  sal_clm(:,:) = reshape(sal_clm_ij, (/nx,ny/) )
1888  endif
1889 
1890 end subroutine get_sal_clm
1891 
1904 subroutine get_sal_clm_ta(sal_clm_ta,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
1907  implicit none
1908 
1909 ! input
1910  integer, intent(in) :: nlat,nlon,mon1,mon2
1911  real, intent(in) :: wei1,wei2
1912 ! output
1913  real, dimension(nlon,nlat), intent(out) :: sal_clm_ta
1914  real, dimension(nlat), intent(out) :: xlats
1915  real, dimension(nlon), intent(out) :: xlons
1916 
1917 !input/output data file names
1918  character (len=6), parameter :: fin_sal_clm='salclm'
1919 
1920 ! local declare
1921  real, dimension(nlon,nlat) :: sal_clm1,sal_clm2
1922 
1923 !
1924 ! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
1925 !
1926  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm1,xlats,xlons,nlat,nlon,mon1)
1927  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm2,xlats,xlons,nlat,nlon,mon2)
1928 !
1929 ! sal_clim at the analysis time
1930 !
1931  sal_clm_ta(:,:) = wei1*sal_clm1(:,:)+wei2*sal_clm2(:,:)
1932  write(*,'(a,2f9.3)') 'sal_clm_ta, min, max : ',minval(sal_clm_ta),maxval(sal_clm_ta)
1933  end subroutine get_sal_clm_ta
1934 
1949 subroutine intp_tile(tf_lalo,dlats_lalo,dlons_lalo,jdim_lalo,idim_lalo, &
1950  tf_tile,xlats_tile,xlons_tile,jdim_tile,idim_tile)
1952  use utils
1953 
1954  implicit none
1955 
1956 ! input/output
1957  real, dimension(idim_lalo,jdim_lalo), intent(in) :: tf_lalo
1958  real, dimension(jdim_lalo), intent(in) :: dlats_lalo
1959  real, dimension(idim_lalo), intent(in) :: dlons_lalo
1960  real, dimension(jdim_tile*idim_tile), intent(in) :: xlats_tile
1961  real, dimension(jdim_tile*idim_tile), intent(in) :: xlons_tile
1962  integer, intent(in) :: jdim_lalo,idim_lalo,jdim_tile,idim_tile
1963  real, dimension(jdim_tile*idim_tile), intent(out) :: tf_tile
1964 
1965 ! local
1966  real, parameter :: deg2rad=3.1415926/180.0
1967  real, dimension(jdim_lalo) :: xlats_lalo
1968  real, dimension(idim_lalo) :: xlons_lalo
1969  real :: wsum
1970  integer :: itile,jtile
1971  integer :: ij
1972  integer :: ilalo,jlalo,ilalop1,jlalop1
1973 
1974  integer, allocatable, dimension(:,:) :: id1,id2,jdc
1975  real, allocatable, dimension(:,:,:) :: agrid,s2c
1976 
1977  print*
1978  print*,'interpolate from lat/lon grids to any one grid with known lat/lon'
1979 
1980  xlats_lalo = dlats_lalo*deg2rad
1981  xlons_lalo = dlons_lalo*deg2rad
1982 
1983  allocate(agrid(idim_tile,jdim_tile,2))
1984  agrid(:,:,1) = reshape(xlons_tile, (/idim_tile,jdim_tile/) )
1985  agrid(:,:,2) = reshape(xlats_tile, (/idim_tile,jdim_tile/) )
1986  agrid = agrid*deg2rad
1987 
1988  allocate(id1(idim_tile,jdim_tile))
1989  allocate(id2(idim_tile,jdim_tile))
1990  allocate(jdc(idim_tile,jdim_tile))
1991  allocate(s2c(idim_tile,jdim_tile,4))
1992 
1993 !----------------------------------------------------------------------
1994 ! compute bilinear weights for each model point from the nearest
1995 ! four lalo points. does not account for mask. that
1996 ! happens later.
1997 !----------------------------------------------------------------------
1998 
1999  call remap_coef( 1, idim_tile, 1, jdim_tile, idim_lalo, jdim_lalo, &
2000  xlons_lalo, xlats_lalo, id1, id2, jdc, s2c, agrid )
2001 
2002  do ij = 1, jdim_tile*idim_tile
2003 
2004  jtile = (ij-1)/idim_tile + 1
2005  itile = mod(ij,idim_tile)
2006  if (itile==0) itile = idim_tile
2007 
2008  ilalo = id1(itile,jtile)
2009  ilalop1 = id2(itile,jtile)
2010  jlalo = jdc(itile,jtile)
2011  jlalop1 = jdc(itile,jtile) + 1
2012 
2013  wsum = s2c(itile,jtile,1) + s2c(itile,jtile,2) + &
2014  s2c(itile,jtile,3) + s2c(itile,jtile,4)
2015 
2016  tf_tile(ij) = ( s2c(itile,jtile,1)*tf_lalo(ilalo,jlalo) + &
2017  s2c(itile,jtile,2)*tf_lalo(ilalop1,jlalo) + &
2018  s2c(itile,jtile,3)*tf_lalo(ilalop1,jlalop1) + &
2019  s2c(itile,jtile,4)*tf_lalo(ilalo,jlalop1) )/wsum
2020  enddo
2021 
2022  deallocate(id1, id2, jdc, s2c)
2023 
2024 end subroutine intp_tile
2025 
2038 subroutine get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
2039  implicit none
2040 
2041 ! input
2042  integer, intent(in) :: iy,im,id,ih
2043 ! output
2044  integer, intent(out) :: mon1,mon2
2045  real, intent(out) :: wei1,wei2
2046 
2047 ! local declare
2048  real :: rjday
2049  integer :: mon,monend,monm,monp,jdow,jdoy,jday
2050  integer :: jda(8)
2051 !
2052 !dayhf : julian day of the middle of each month
2053 !
2054  real, dimension(13) :: dayhf
2055  data dayhf/15.5,45.0,74.5,105.0,135.5,166.0,196.5,227.5,258.0,288.5,319.0,349.5,380.5/
2056 
2057 ! 15, 44, 73.5, 104, 134.5, 165, 195.5, 226.5, 257, 287.5, 318.5, 349 ! from
2058 ! woa05
2059 
2060  jda=0
2061  jda(1)=iy
2062  jda(2)=im
2063  jda(3)=id
2064  jda(5)=ih
2065  jdow = 0
2066  jdoy = 0
2067  jday = 0
2068  call w3doxdat(jda,jdow,jdoy,jday)
2069  rjday=jdoy+jda(5)/24.
2070  if(rjday.lt.dayhf(1)) rjday=rjday+365.
2071 
2072  monend = 12
2073  do mon = 1, monend
2074  monm = mon
2075  monp = mon + 1
2076  if( rjday >= dayhf(monm) .and. rjday < dayhf(monp) ) then
2077  mon1 = monm
2078  mon2 = monp
2079  go to 10
2080  endif
2081  enddo
2082 
2083  print *,'FATAL ERROR in get_tim_wei, wrong rjday',rjday
2084  call abort
2085  10 continue
2086 
2087  wei1 = (dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
2088  wei2 = (rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))
2089 
2090  if( mon2 == 13 ) mon2=1
2091 
2092  write(*,'(a,2i4,3f9.3)') 'mon1,mon2,rjday,wei1,wei2=',mon1,mon2,rjday,wei1,wei2
2093 
2094  end subroutine get_tim_wei
2095 
2105  real function tfreez(salinity)
2107  implicit none
2108 
2109  REAL salinity,sal
2110  REAL a1, a2, a3
2111  parameter(a1 = -0.0575)
2112  parameter(a2 = 1.710523e-3)
2113  parameter(a3 = -2.154996e-4)
2114 
2115  IF (salinity .LT. 0.) THEN
2116  sal = 0.
2117  ELSE
2118  sal = salinity
2119  ENDIF
2120  tfreez = sal*(a1+a2*sqrt(sal)+a3*sal)
2121 
2122  return
2123  end
integer function num_parthds()
Return the number of omp threads.
Definition: num_parthds.f90:10
subroutine, public read_data(LSOIL, LENSFC, DO_NSST, DO_SNO_INC_JEDI, DO_SOI_INC_JEDI, INC_FILE, IS_NOAHMP, TSFFCS, SMCFCS, SWEFCS, STCFCS, TG3FCS, ZORFCS, CVFCS, CVBFCS, CVTFCS, ALBFCS, VEGFCS, SLIFCS, CNPFCS, F10M, VETFCS, SOTFCS, ALFFCS, USTAR, FMM, FHH, SIHFCS, SICFCS, SITFCS, TPRCP, SRFLAG, SNDFCS, VMNFCS, VMXFCS, SLCFCS, STCINC, SLCINC, SLPFCS, ABSFCS, T2M, Q2M, SLMASK, ZSOIL, NSST)
Read the first guess surface records and nsst records (if selected) for a single cubed-sphere tile...
subroutine, public get_tf_clm_dim(file_sst, mlat_sst, mlon_sst)
Get the i/j dimensions of RTG SST climatology file.
subroutine, public read_tf_clim_grb(file_sst, sst, rlats_sst, rlons_sst, mlat_sst, mlon_sst, mon)
Read a GRIB1 sst climatological analysis file.
integer, public idim_gaus
&#39;i&#39; dimension of GSI gaussian grid.
subroutine get_tf_clm(xlats_ij, xlons_ij, ny, nx, iy, im, id, ih, tf_clm, tf_trd)
Get the sst climatology at the valid time and on the target grid.
Definition: cycle.f90:1728
subroutine, public read_salclm_gfs_nc(filename, sal, xlats, xlons, nlat, nlon, itime)
Read the woa05 salinity monthly climatology file.
subroutine get_tf_clm_ta(tf_clm_ta, tf_clm_trend, xlats, xlons, nlat, nlon, mon1, mon2, wei1, wei2)
Get the reference temperature/sst climatology and its trend at analysis time.
Definition: cycle.f90:1797
subroutine, public remap_coef(is, ie, js, je, im, jm, lon, lat, id1, id2, jdc, s2c, agrid)
Generate the weights and index of the grids used in the bilinear interpolation.
Definition: utils.F90:40
integer, dimension(:,:), allocatable, public slmsk_gaus
GSI land mask on the gaussian grid.
subroutine dtzm_point(XT, XZ, DT_COOL, ZC, Z1, Z2, DTZM)
Compute the vertical mean of the NSST t-profile.
Definition: cycle.f90:1519
subroutine get_sal_clm_ta(sal_clm_ta, xlats, xlons, nlat, nlon, mon1, mon2, wei1, wei2)
Get climatological salinity at the analysis time.
Definition: cycle.f90:1905
Holds machine dependent constants for global_cycle.
Definition: machine.f90:7
subroutine sfcdrv(LUGB, IDIM, JDIM, LENSFC, LSOIL, DELTSFC, IY, IM, ID, IH, FH, IALB, USE_UFO, DO_NSST, DO_SFCCYCLE, DO_LNDINC, FRAC_GRID, ZSEA1, ZSEA2, ISOT, IVEGSRC, MYRANK)
Driver routine for updating the surface fields.
Definition: cycle.f90:321
subroutine tf_thaw_set(tf_ij, mask_ij, itile, jtile, tice, tclm, tf_thaw, nx, ny, nset_thaw_s, nset_thaw_i, nset_thaw_c)
Set the background reference temperature (tf) for points where the ice has just melted.
Definition: cycle.f90:1592
subroutine, public get_dim_nc(filename, nlat, nlon)
Get the i/j dimensions of the data from a NetCDF file.
subroutine, public read_gsi_data(GSI_FILE, FILE_TYPE, LSOIL)
Read increment file from the GSI containing either the foundation temperature increments and mask...
subroutine intp_tile(tf_lalo, dlats_lalo, dlons_lalo, jdim_lalo, idim_lalo, tf_tile, xlats_tile, xlons_tile, jdim_tile, idim_tile)
Interpolate lon/lat grid data to the fv3 native grid (tf_lalo => tf_tile).
Definition: cycle.f90:1951
subroutine climo_trend(LATITUDE, MON, DAY, DELTSFC, DTREF)
If the tile point is an isolated water point that has no corresponding gsi water point, then tref is updated using the rtg sst climo trend.
Definition: cycle.f90:1366
real, dimension(:,:), allocatable, public dtref_gaus
GSI foundation temperature increment on the gaussian grid.
subroutine, public write_data(lensfc, idim, jdim, lsoil, do_nsst, inc_file, nsst, slifcs, tsffcs, vegfcs, swefcs, tg3fcs, zorfcs, albfcs, alffcs, cnpfcs, f10m, t2m, q2m, vetfcs, sotfcs, ustar, fmm, fhh, sicfcs, sihfcs, sitfcs, tprcp, srflag, swdfcs, vmnfcs, vmxfcs, slpfcs, absfcs, slcfcs, smcfcs, stcfcs, stcinc, slcinc)
Update surface records - and nsst records if selected - on a single cubed-sphere tile to a pre-existi...
Module containing utility routines.
Definition: utils.F90:7
This module contains routines that read and write data.
program sfc_drv
Stand alone surface/NSST cycle driver for the cubed-sphere grid.
Definition: cycle.f90:111
subroutine adjust_nsst(RLA, RLO, SLMSK_TILE, SLMSK_FG_TILE, SKINT_TILE, SICET_TILE, sice_tile, sice_fg_tile, SOILT_TILE, NSST, LENSFC, LSOIL, IDIM, JDIM, ZSEA1, ZSEA2, tf_clm_tile, tf_trd_tile, sal_clm_tile, LANDFRAC, FRAC_GRID)
Read in gsi file with the updated reference temperature increments (on the gaussian grid)...
Definition: cycle.f90:942
real function tfreez(salinity)
Compute the freezing point of water as a function of salinity.
Definition: cycle.f90:2106
integer, public jdim_gaus
&#39;j&#39; dimension of GSI gaussian grid.
subroutine get_tim_wei(iy, im, id, ih, mon1, mon2, wei1, wei2)
For a given date, determine the bounding months and the linear time interpolation weights...
Definition: cycle.f90:2039
subroutine get_sal_clm(xlats_ij, xlons_ij, ny, nx, iy, im, id, ih, sal_clm)
Get salinity climatology at the valid time on the target grid.
Definition: cycle.f90:1845
subroutine nsst_water_reset(nsst, ij, tf_thaw)
If the first guess was sea ice, but the analysis is open water, reset all nsst variables.
Definition: cycle.f90:1682
subroutine, public read_lat_lon_orog(RLA, RLO, OROG, OROG_UF, TILE_NUM, IDIM, JDIM, IJDIM, LANDFRAC)
Read latitude and longitude for the cubed-sphere tile from the &#39;grid&#39; file.