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