MODULE AtmosModel
  USE Constants, ONLY: &
       r8,             & 
       r4,             &
       i8,             &
       pi
  USE Utils, ONLY :julday,lati
  USE Options, ONLY :lon_site, lat_site,icn_data,DELTAOUT
  !PK USE ReadNETCDF    , ONLY : RunReadNETCDF,ncid_prec,ncid_temp,ncid_rhum,ncid_wind,ncid_swrd,ncid_lwrd,ncid_umes,ncid_pres

  IMPLICIT NONE
  PRIVATE

  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gu	  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gv	  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gt	  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gq	  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gps	  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: dcupr_sib  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: dlspr_sib  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: radvbc_sib (:  ,:)!calc em funcao de swdown
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: radnbc_sib (:  ,:)!calc em funcao de swdown
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: radvdc_sib (:  ,:)!calc em funcao de swdown
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: radndc_sib (:  ,:)!calc em funcao de swdown
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: swdown(:,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: dlwbot_sib (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: cosz_sib   (:  ,:)!calc em funcao de latitude
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: ws	  (:,:,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: psb_g	  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: bps_g	  (:  ,:)
  REAL(KIND=r8), PUBLIC            :: sigki
  REAL(KIND=r8), PUBLIC            :: delsig
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: mskant8    (:  ,:)  
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: tseam	  (:  ,:)    
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: tsea	  (:  ,:)    
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: slrad	  (:  ,:)    
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: tsurf	  (:  ,:)    
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: qsurf	  (:  ,:)    
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: zorl	  (:  ,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: tmtx	  (:,:,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: qmtx	  (:,:,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: umtx	  (:,:,:)
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: slhf	  (:  ,:)    
  REAL(KIND=r8), PUBLIC, ALLOCATABLE :: sshf	  (:  ,:)    

  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: latic(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: longc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: t02mc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: q02mc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: uvesc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: vvesc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: pslcc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: precc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: ocisc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: slhfc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: sshfc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: olisc(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: latip(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: longp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: t02mp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: q02mp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: uvesp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: vvesp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: pslcp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: precp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: ocisp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: slhfp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: sshfp(:  ,:)
  REAL(KIND=r4), PUBLIC, ALLOCATABLE :: olisp(:  ,:)
  INTEGER :: monday(12)
  INTEGER           , PUBLIC           :: monl(12)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
  REAL (KIND=r8), PARAMETER   :: pie  = 	     3.1415926e0_r8! constant pi=3.1415926e0   
  REAL (KIND=r8), PARAMETER   :: pai12  =	       pie/12.0_r8 ! pi/12
  REAL (KIND=r8), PARAMETER   :: pai   =        3.14159265358979_r8! constant pi=3.1415926e0
  LOGICAL :: OPENFILE=.TRUE.
  INTEGER :: ibMax
  INTEGER :: jbMax
  INTEGER :: nband
  INTEGER :: step

  PUBLIC :: InitAtmosModel
  PUBLIC :: RunAtmosModel2D
  PUBLIC :: FinalizeAtmosModel
  PUBLIC :: RunAtmosModel1D
  
CONTAINS

  SUBROUTINE InitAtmosModel(ibMax_in,jbMax_in,nband_in)
    IMPLICIT NONE
    INTEGER, INTENT(IN  ) :: ibMax_in
    INTEGER, INTENT(IN  ) :: jbMax_in
    INTEGER, INTENT(IN  ) :: nband_in
    ibMax = ibMax_in
    jbMax = jbMax_in
    nband = nband_in
    ALLOCATE(gu         (1:ibMax    ,jbMax)  )
    ALLOCATE(gv         (1:ibMax    ,jbMax)  )
    ALLOCATE(gt         (1:ibMax    ,jbMax)  )
    ALLOCATE(gq         (1:ibMax    ,jbMax)  )
    ALLOCATE(gps        (1:ibMax    ,jbMax)  )
    ALLOCATE(dcupr_sib  (1:ibMax    ,jbMax)  )
    ALLOCATE(dlspr_sib  (1:ibMax    ,jbMax)  )
    ALLOCATE(radvbc_sib (1:ibMax    ,jbMax)  )
    ALLOCATE(radnbc_sib (1:ibMax    ,jbMax)  )
    ALLOCATE(radvdc_sib (1:ibMax    ,jbMax)  )
    ALLOCATE(radndc_sib (1:ibMax    ,jbMax)  )
    ALLOCATE(dlwbot_sib (1:ibMax    ,jbMax)  )
    ALLOCATE(swdown     (1:ibMax    ,jbMax)  )
    ALLOCATE(cosz_sib   (1:ibMax    ,jbMax)  )
    ALLOCATE(ws         (1:ibMax,1:3,jbMax)  )
    ALLOCATE(psb_g      (1:ibMax    ,jbMax)  )
    ALLOCATE(bps_g	(1:ibMax    ,jbMax)  )
    ALLOCATE(mskant8    (1:ibMax    ,jbMax)  )
    ALLOCATE(tseam      (1:ibMax    ,jbMax)  )
    ALLOCATE(tsea       (1:ibMax    ,jbMax)  )
    ALLOCATE(slrad      (1:ibMax    ,jbMax)  )
    ALLOCATE( slhf     (1:ibMax    ,jbMax)  )
    ALLOCATE( sshf     (1:ibMax    ,jbMax)  )
    ALLOCATE(tsurf      (1:ibMax    ,jbMax)  )
    ALLOCATE(qsurf      (1:ibMax    ,jbMax)  )
    ALLOCATE(zorl       (1:ibMax    ,jbMax)  )
    ALLOCATE(tmtx       (1:ibMax,1:3,jbMax)  )
    ALLOCATE(qmtx	     (1:ibMax,1:3,jbMax)  )
    ALLOCATE(umtx	     (1:ibMax,1:4,jbMax)  )

    ALLOCATE(latic(1:ibMax    ,jbMax)  )
    ALLOCATE(longc(1:ibMax    ,jbMax)  )
    ALLOCATE(t02mc(1:ibMax    ,jbMax)  )
    ALLOCATE(q02mc(1:ibMax    ,jbMax)  )
    ALLOCATE(uvesc(1:ibMax    ,jbMax)  )
    ALLOCATE(vvesc(1:ibMax    ,jbMax)  )
    ALLOCATE(pslcc(1:ibMax    ,jbMax)  )
    ALLOCATE(precc(1:ibMax    ,jbMax)  )
    ALLOCATE(ocisc(1:ibMax    ,jbMax)  )
    ALLOCATE(slhfc(1:ibMax    ,jbMax)  )
    ALLOCATE(sshfc(1:ibMax    ,jbMax)  )
    ALLOCATE(olisc(1:ibMax    ,jbMax)  )
    ALLOCATE(latip(1:ibMax    ,jbMax)  )
    ALLOCATE(longp(1:ibMax    ,jbMax)  )
    ALLOCATE(t02mp(1:ibMax    ,jbMax)  )
    ALLOCATE(q02mp(1:ibMax    ,jbMax)  )
    ALLOCATE(uvesp(1:ibMax    ,jbMax)  )
    ALLOCATE(vvesp(1:ibMax    ,jbMax)  )
    ALLOCATE(pslcp(1:ibMax    ,jbMax)  )
    ALLOCATE(precp(1:ibMax    ,jbMax)  )
    ALLOCATE(ocisp(1:ibMax    ,jbMax)  )
    ALLOCATE(slhfp(1:ibMax    ,jbMax)  )
    ALLOCATE(sshfp(1:ibMax    ,jbMax)  )
    ALLOCATE(olisp(1:ibMax    ,jbMax)  )

    CALL InitRadtim(monl)
  END SUBROUTINE InitAtmosModel

  SUBROUTINE FinalizeAtmosModel()
    IMPLICIT NONE
    DEALLOCATE(gu        )
    DEALLOCATE(gv        )
    DEALLOCATE(gt        )
    DEALLOCATE(gq        )
    DEALLOCATE(gps       )
    DEALLOCATE(dcupr_sib )
    DEALLOCATE(dlspr_sib )
    DEALLOCATE(radvbc_sib)
    DEALLOCATE(radnbc_sib)
    DEALLOCATE(radvdc_sib)
    DEALLOCATE(radndc_sib)
    DEALLOCATE(dlwbot_sib)
    DEALLOCATE(    swdown)
    DEALLOCATE(cosz_sib  )
    DEALLOCATE(ws        )
    DEALLOCATE(psb_g     )
    DEALLOCATE(bps_g     )
    DEALLOCATE(mskant8   )
    DEALLOCATE(tseam     )
    DEALLOCATE(tsea      )
    DEALLOCATE(slrad     )
    DEALLOCATE( slhf    )
    DEALLOCATE( sshf    )
    DEALLOCATE(tsurf     )
    DEALLOCATE(qsurf     )
    DEALLOCATE(zorl      )
    DEALLOCATE(tmtx      )
    DEALLOCATE(qmtx      )
    DEALLOCATE(umtx      )
    DEALLOCATE(latic     )
    DEALLOCATE(longc     )
    DEALLOCATE(t02mc     )
    DEALLOCATE(q02mc     )
    DEALLOCATE(uvesc     )
    DEALLOCATE(vvesc     )
    DEALLOCATE(pslcc     )
    DEALLOCATE(precc     )
    DEALLOCATE(ocisc     )
    DEALLOCATE(slhfc     )
    DEALLOCATE(sshfc     )
    DEALLOCATE(olisc     )
    DEALLOCATE(latip     )
    DEALLOCATE(longp     )
    DEALLOCATE(t02mp     )
    DEALLOCATE(q02mp     )
    DEALLOCATE(uvesp     )
    DEALLOCATE(vvesp     )
    DEALLOCATE(pslcp     )
    DEALLOCATE(precp     )
    DEALLOCATE(ocisp     )
    DEALLOCATE(slhfp     )
    DEALLOCATE(sshfp     )
    DEALLOCATE(olisp     )

  END SUBROUTINE FinalizeAtmosModel

  SUBROUTINE RunAtmosModel1D(jdt,tod,nband,dt,idate,idatec,idatep,DELTAIN,start)
    IMPLICIT NONE
    REAL(KIND=r8), INTENT(IN   ) :: tod
    INTEGER       , INTENT(IN   ) :: jdt
    INTEGER       , INTENT(IN   ) :: nband  ! INTENT(IN    ) :: nband
    REAL(KIND=r8), INTENT(IN   ) :: dt     ! INTENT(IN    ) :: dt 
    INTEGER       , INTENT(IN   ) :: idate (4) ! INTENT(IN    ) :: idatec(:)
    INTEGER       , INTENT(IN   ) :: idatec(4) ! INTENT(IN    ) :: idatec(:)
    INTEGER       , INTENT(IN   ) :: idatep(4)
    REAL(KIND=r8), INTENT(IN   ) :: DELTAIN
    LOGICAL       , INTENT(IN   ) :: start
    INTEGER :: nCols
    INTEGER        :: jb
    INTEGER        :: ib
    INTEGER        :: i
    INTEGER        :: ios
    REAL(KIND=r8) :: rtime
    REAL(KIND=r8) :: time
    REAL(KIND=r8) :: time2
    REAL(KIND=r8) :: jday
    REAL(KIND=r8) :: orbit
    REAL(KIND=r8) :: angle
    REAL(KIND=r8) :: xdecl
    REAL(KIND=r8) :: frac
    REAL(KIND=r8) :: sw
    REAL(KIND=r8) :: year
    REAL(KIND=r8) :: decmax
    REAL(KIND=r8) :: sols
    REAL(KIND=r8) :: season
    REAL(KIND=r8) :: day
    REAL(KIND=r8) :: dec
    REAL(KIND=r8) :: rfd
    REAL(KIND=r8) :: sind
    REAL(KIND=r8) :: cosd
    REAL(KIND=r8) :: hac
    REAL(KIND=r8) :: h
    REAL(KIND=r8) :: sr
    REAL(KIND=r8) :: ss
    REAL(KIND=r8) :: dawn 
    REAL(KIND=r8) :: dusk  
    REAL(KIND=r8) :: coshr 
    REAL(KIND=r8) :: sunang
    REAL(KIND=r8) :: trans
    REAL(KIND=r8) :: btime 
    REAL(KIND=r8) :: ctime
    REAL(KIND=r8) :: atime 
    REAL(KIND=r8) :: fdiffuse
    REAL(KIND=r8) :: xlati (jbMax)
    REAL(KIND=r8) :: coszen(ibMax,jbMax)
    REAL(KIND=r8) :: cloud(ibMax,jbMax)
    REAL(KIND=r8) :: solad(ibMax,nband)
    REAL(KIND=r8) :: solai(ibMax,nband)
    REAL(KIND=r8) :: var(ibMax,jbMax)
    CHARACTER(LEN=255) :: namec
    CHARACTER(LEN=255) :: namep
    INTEGER            :: open_status
    INTEGER            :: reclen
    INTEGER            :: irec
    INTEGER            :: j
    INTEGER, PARAMETER   :: input_unit=200
    REAL(KIND=r8) , PARAMETER    :: yrl=365.2500_r8
    INTEGER, PARAMETER   :: ok = 0
    REAL(KIND=r8)   , PARAMETER :: alon   =  0.0_r8
    REAL (KIND=r8), PARAMETER   :: pie  =              3.1415926e0_r8! constant pi=3.1415926e0   
    REAL (KIND=r8), PARAMETER   :: pai12  =              pie/12.0_r8 ! pi/12

    LOGICAL :: here
    CHARACTER(LEN=255) :: PREFIX='GCM2IBIS'
    CHARACTER(LEN=255) :: SUFFIX='.dat'
    INTEGER           :: iyeard 
    INTEGER           :: imonth
    INTEGER           :: juliand
    REAL(KIND=r8)    :: frachours
    REAL(KIND=r8)    :: hours
    REAL(KIND=r8)    :: xua
    REAL(KIND=r8)    :: xta
    REAL(KIND=r8)    :: xprec
    REAL(KIND=r8)    :: xsin
    REAL(KIND=r8)    :: xlin
    REAL(KIND=r8)    :: xur
    REAL(KIND=r8)    :: xes
    REAL(KIND=r8)    :: xe
    REAL(KIND=r8)    :: xqa
    REAL(KIND=r8)    :: xpslc
    REAL(KIND=r8)    :: xlh
    REAL(KIND=r8)    :: xsh
    REAL(KIND=r8)    :: sdelt
    REAL(KIND=r8)    :: ratio
    REAL(KIND=r8)    :: etime
    REAL(KIND=r8)    :: xday
    REAL(KIND=r8)    :: sindel
    REAL(KIND=r8)    :: cosdel
    REAL(KIND=r8)    :: fimxi 
    REAL(KIND=r8)    :: frh
    REAL(KIND=r8)    :: cos2  (ibMax)
    REAL(KIND=r8)    :: zenith1 (ibMax)
    REAL(KIND=r8)    :: zenith2 (ibMax)
    REAL(KIND=r8)    :: lonrad(ibMax)
    REAL(KIND=r8)    :: swradtot    (ibMax,jbMax) 
    REAL(KIND=r8)    :: diffswradtot(ibMax,jbMax)

    open_status=ok
    nCols=ibMax
    time=tod
    IF(MOD(tod,DELTAIN) == 0.0_r8 .OR. jdt == 1)THEN
       step=0
       DO i=1,ibMax
          WRITE(namec,'(A12,i3.3)')'Clim_input.S',i    
          INQUIRE (FILE='./inputdata/'//TRIM(namec),EXIST=here)
          IF (.NOT.here) THEN
             WRITE(0,*)'FILE ','./inputdata/'//TRIM(namec), 'Not Exist';STOP
          ELSE
             IF(OPENFILE)OPEN(input_unit+i,FILE='./inputdata/'//TRIM(namec),ACCESS='SEQUENTIAL',&
                FORM='FORMATTED',STATUS='OLD',ACTION='READ', &
                IOSTAT=open_status)
             IF( open_status == ok ) THEN
                IF(TRIM(icn_data) == 'ERA5')THEN
                   READ(input_unit+i,'(3I8.8,9F8.3)') iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xlh,xsh
                ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
                   READ(input_unit+i,*)               iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xpslc,xlh,xsh
                ELSE
                   PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                   STOP       
                END IF
                DO j=1,jbMax
                
                   latic(i,j) = lat_site(i)!degree
                   longc(i,j) = lon_site(i)!degree
                   t02mc(i,j) = xta+273.16_r8!K
                   IF(TRIM(icn_data) == 'ERA5')THEN
                      xpslc      =( 1.15_r8*287.16*(xta+273.16_r8)  )/100.0_r8   ! (rho*R*T)  hPa
                      xes=6.112*exp(17.67*xta/(243.5+xta))!hPa
                      xe=(xur/100.0_r8)*xes     !hPa
                      xqa  =xe*622.0/(xpslc-xe) !g/kg
                      q02mc(i,j) = xqa/1000.0_r8!kg/kg
                   ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
                      q02mc(i,j) = xur!kg/kg
                   ELSE 
                      PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                      STOP       
                   END IF

                   uvesc(i,j) = xua!m/s
                   vvesc(i,j) = xua!m/s
                   pslcc(i,j) = xpslc!hPa
                   IF(TRIM(icn_data) == 'ERA5')THEN
                       precc(i,j) = xprec*86400/DELTAIN
                       ocisc(i,j) = xsin!w/m2
                       olisc(i,j) = xlin!w/m2
                       slhfc(i,j) = xlh
                       sshfc(i,j) = xsh
                   ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
                       precc(i,j) = xprec 
                       ocisc(i,j) = xsin
                       olisc(i,j) = xlin
                       slhfc(i,j) = xlh
                       sshfc(i,j) = xsh
                   ELSE 
                      PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                      STOP       
                   END IF
                END DO
                   
                IF(.FALSE.)THEN
                   WRITE(0,*)'latic',MAXVAL(latic),MINVAL(latic)
                   WRITE(0,*)'longc',MAXVAL(longc),MINVAL(longc)
                   WRITE(0,*)'t02mc',MAXVAL(t02mc),MINVAL(t02mc)
                   WRITE(0,*)'q02mc',MAXVAL(q02mc),MINVAL(q02mc)
                   WRITE(0,*)'uvesc',MAXVAL(uvesc),MINVAL(uvesc)
                   WRITE(0,*)'vvesc',MAXVAL(vvesc),MINVAL(vvesc)
                   WRITE(0,*)'pslcc',MAXVAL(pslcc),MINVAL(pslcc)
                   WRITE(0,*)'precc',MAXVAL(precc),MINVAL(precc)
                   WRITE(0,*)'ocisc',MAXVAL(ocisc),MINVAL(ocisc)
                   WRITE(0,*)'olisc',MAXVAL(olisc),MINVAL(olisc)
                   WRITE(0,*)'slhfc',MAXVAL(slhfc),MINVAL(slhfc)
                   WRITE(0,*)'sshfc',MAXVAL(sshfc),MINVAL(sshfc)
                END IF
	     ELSE
                WRITE(0,*)"Unable to OPEN file";STOP
             END IF
          END IF
          IF(start)REWIND(input_unit+i)
       END DO 

       DO i=1,ibMax
          WRITE(namep,'(A13,i3.3)')'Clima_input.S',i    
       
          INQUIRE (FILE='./inputdata/'//TRIM(namep),EXIST=here)
          IF (.NOT.here) THEN
             WRITE(0,*)'FILE ','./inputdata/'//TRIM(namep), 'Not Exist';STOP
          ELSE
             INQUIRE(IOLENGTH=reclen)latic
             IF(OPENFILE)OPEN(input_unit+ibMax+i,FILE='./inputdata/'//TRIM(namep),ACCESS='SEQUENTIAL',&
                FORM='FORMATTED',STATUS='OLD',ACTION='READ', &
                IOSTAT=open_status)
             IF( open_status == ok ) THEN
                IF(TRIM(icn_data) == 'ERA5')THEN
                   READ(input_unit+ibMax+i,'(3I8.8,9F8.3)') iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xlh,xsh
                ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
                   IF(OPENFILE)READ(input_unit+ibMax+i,*)               iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xpslc,xlh,xsh
                   READ(input_unit+ibMax+i,*)               iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xpslc,xlh,xsh
                ELSE
                   PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                   STOP       
                END IF
                DO j=1,jbMax		
                   latip(i,j) = lat_site(i)!degree
                   longp(i,j) = lon_site(i)!degree
                   t02mp(i,j) = xta+273.16_r8!K
                  IF(TRIM(icn_data) == 'ERA5')THEN
                      xpslc      = (1.15_r8*287.16*(xta+273.16_r8))/100.0_r8     ! (rho*R*T)  hPa
                      xes=6.112*exp(17.67*xta/(243.5+xta))!hPa
                      xe=(xur/100.0_r8)*xes     !hPa
                      xqa  =xe*622.0/(xpslc-xe) !g/kg
                      q02mp(i,j) = xqa/1000.0_r8!kg/kg
                   ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
                      q02mc(i,j) = xur!kg/kg
                   ELSE  
                      PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                      STOP       
                   END IF
                   uvesp(i,j) = xua!m/s
                   vvesp(i,j) = xua!m/s
                   pslcp(i,j) = xpslc!hPaxpslc!hPa
                   IF(TRIM(icn_data) == 'ERA5')THEN
                       precp(i,j) = xprec*86400/DELTAIN
                       ocisp(i,j) = xsin!w/m2
                       olisp(i,j) = xlin!w/m2
                       slhfp(i,j) = xlh
                       sshfp(i,j) = xsh
                   ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
                       precp(i,j) = xprec  
                       ocisp(i,j) = xsin
                       olisp(i,j) = xlin
                       slhfp(i,j) = xlh 
                       sshfp(i,j) = xsh 
                   ELSE 
                      PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                      STOP       
                   END IF
		END DO
                IF(.FALSE.)THEN
                   WRITE(0,*)'latip',MAXVAL(latip),MINVAL(latip)
                   WRITE(0,*)'longp',MAXVAL(longp),MINVAL(longp)
                   WRITE(0,*)'t02mp',MAXVAL(t02mp),MINVAL(t02mp)
                   WRITE(0,*)'q02mp',MAXVAL(q02mp),MINVAL(q02mp)
                   WRITE(0,*)'uvesp',MAXVAL(uvesp),MINVAL(uvesp)
                   WRITE(0,*)'vvesp',MAXVAL(vvesp),MINVAL(vvesp)
                   WRITE(0,*)'pslcp',MAXVAL(pslcp),MINVAL(pslcp)
                   WRITE(0,*)'precp',MAXVAL(precp),MINVAL(precp)
                   WRITE(0,*)'ocisp',MAXVAL(ocisp),MINVAL(ocisp)
                   WRITE(0,*)'olisp',MAXVAL(olisp),MINVAL(olisp)
                   WRITE(0,*)'slhfp',MAXVAL(slhfp),MINVAL(slhfp)
                   WRITE(0,*)'sshfp',MAXVAL(sshfp),MINVAL(sshfp)
                END IF

             ELSE
                WRITE(0,*)"Unable to OPEN file";STOP
             END IF
          END IF
          IF(start)REWIND(input_unit+ibMax+i)
       END DO     
       WRITE(0,*)'READ ',TRIM(namec),' and ',TRIM(namep),iyeard,imonth
    END IF
    OPENFILE=.FALSE.
    DO jb=1,jbMax
        !pi --- 180
	! y      x
	!
	! y =  X*pi/180 !radianos
	! 
        lati   (jbMax+1-jb)    = (latic(1,jb)+90.0_r8)*pi/180.0_r8
        xlati  (jbMax+1-jb)    = (latic(1,jb))*pi/180.0_r8
       DO ib=1,nCols
          lonrad(ib)                = longc(ib,1)
          gu        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*uvesc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*uvesp(ib,jb)
          gv        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*vvesc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*vvesp(ib,jb)
          gt        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          gq        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*q02mp(ib,jb)
          gps       (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*pslcc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*pslcp(ib,jb)
          var       (ib,        jb) = ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*precp(ib,jb)

          IF(TRIM(icn_data) == 'ERA5')THEN
              dcupr_sib (ib,jbMax+1-jb) = (var(ib,jb)/86400.0_r8)*1.0e+3_r8
             !          PRINT*,  var       (ib,        jb), ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb), ((step*dt)/(DELTAIN))*precp(ib,jb)
              dlspr_sib (ib,jbMax+1-jb) = 0.0_r8
          ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
              dcupr_sib (ib,jbMax+1-jb) = var(ib,jb)*1.0e+3_r8!/ convert m to mm
!             PRINT*,  var       (ib,        jb), ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb), ((step*dt)/(DELTAIN))*precp(ib,jb)
              dlspr_sib (ib,jbMax+1-jb) = 0.0_r8
          ELSE
             PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
             STOP       
           END IF
          swdown    (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*ocisc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*ocisp(ib,jb)
          dlwbot_sib(ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*olisp(ib,jb)
          sigki                     = 1.0_r8
          delsig                    = 0.010000_r8
          tseam     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          tsea      (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          slrad     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*olisp(ib,jb)
          slhf      (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*slhfc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*slhfp(ib,jb)
          sshf      (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*sshf(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*sshfp(ib,jb)
    ! zorl........zorl (i)= 100.0 *zgrav*speedm(i)*rhi(i)
    !             zgrav =0.032 /grav

          zorl      (ib,jbMax+1-jb) = 100.0_r8*(0.32_r8/9.8_r8)*(sqrt( gu(ib,jbMax+1-jb)**2 + gv(ib,jbMax+1-jb)**2))*0.01
          tsurf     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          qsurf     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) &
	                            + ((step*dt)/(DELTAIN))*q02mp(ib,jb)
          tmtx      (ib,:,jbMax+1-jb) = 0.010000_r8
          qmtx      (ib,:,jbMax+1-jb) = 0.010000_r8
          umtx      (ib,:,jbMax+1-jb) = 0.010000_r8


                IF(.TRUE.)THEN
                   WRITE(0,*)'latip',MAXVAL(latip),MINVAL(latip)
                   WRITE(0,*)'longp',MAXVAL(longp),MINVAL(longp)
                   WRITE(0,*)'t02mp',MAXVAL(gt),MINVAL(gt)
                   WRITE(0,*)'q02mp',MAXVAL(gq),MINVAL(gq)
                   WRITE(0,*)'uvesp',MAXVAL(gu),MINVAL(gu)
                   WRITE(0,*)'vvesp',MAXVAL(gv),MINVAL(gv)
                   WRITE(0,*)'pslcp',MAXVAL(gps),MINVAL(gps)
                   WRITE(0,*)'precp',MAXVAL(dcupr_sib),MINVAL(dcupr_sib)
                   WRITE(0,*)'ocisp',MAXVAL(swdown),MINVAL(swdown)
                   WRITE(0,*)'olisp',MAXVAL(dlwbot_sib),MINVAL(dlwbot_sib)
                   WRITE(0,*)'slhfp',MAXVAL(slhfp),MINVAL(slhfp)
                   WRITE(0,*)'sshfp',MAXVAL(sshfp),MINVAL(sshfp)
                END IF

       END DO
    END DO

    step=step+1
    

    !
    !     mon is the month used for vegetation data input
    !
    DO jb=1,jbMax
       !DO ib=1,ibMax,1
       !   month(ib,jb)=idatec(2)
       !   IF((((lati(jb)*180.0_r8)/3.1415926e0_r8)-90.0_r8)  > 0.0_r8 ) THEN
       !      month(ib,jb)  =  month(ib,jb) + 6
       !      IF(month(ib,jb).GE.13) month(ib,jb) = month(ib,jb)-12
       !   END IF
       !END DO
       !
       !     computation of astronomical parameters
       !     sdelt ;solar inclination
       !     etime ;correction factor to local time
       !     ratio ;factor relating to the distance between the earth and the sun
       !
       CALL radtim(idatec,sdelt ,ratio ,etime ,tod   ,xday  ,yrl)
       sw=1370.0_r8*ratio

       sindel = SIN(sdelt)
       cosdel = COS(sdelt)
       fimxi  = 24.0e0_r8 /360.0_r8
       ctime  = alon/15.0e0_r8
       cos2   = 0.0e0_r8
       frh    = ( MOD(tod+0.03125_r8,3600.0_r8)-0.03125_r8)/3600.0_r8
       DO ib=1,ibMax
          zenith1(ib)  = sindel*COS(lati(jb))
       ENDDO
       DO ib=1,ibMax
          btime       = fimxi*lonrad(ib)+ctime
          atime       = etime+pai12*(12.0_r8-idatec(1)-frh-btime)
          zenith2 (ib) = cosdel*SIN(lati(jb))*COS(atime)
          coszen  (ib,jb) =  MAX (0.0_r8,zenith1(ib) + zenith2(ib))
	  
	  sunang = MAX( 0.01_r8, coszen  (ib,jb) )
	  
          cloud(ib,jb) = (sw * sunang - swdown(ib,jb)) / (sw * sunang)
          cloud(ib,jb) = MAX(cloud(ib,jb),0.0_r8)
          cloud(ib,jb) = MIN(cloud(ib,jb),1.0_r8)
          !cloud(ib,jb) = MAX(0.58_r8     ,cloud(ib,jb))
       END DO

    END DO
    !
    ! do for each waveband
    !
    DO ib = 1, nband
       DO jb=1,jbMax
          DO i=1,nCols
             ! calculate the solar transmission through the atmosphere
             ! using simple linear function of tranmission and cloud cover
             !
             ! note that the 'cloud cover' data is typically obtained from
             ! sunshine hours -- not direct cloud observations
             !
             ! where, cloud cover = 1 - sunshine fraction 
             !
             ! different authors present different values for the slope and 
             ! intercept terms of this equation
             !
             ! Friend, A: Parameterization of a global daily weather generator for
             ! terrestrial ecosystem and biogeochemical modelling, Ecological 
             ! Modelling
             !
             ! Spitters et al., 1986: Separating the diffuse and direct component
             ! of global radiation and its implications for modeling canopy
             ! photosynthesis, Part I: Components of incoming radiation,
             ! Agricultural and Forest Meteorology, 38, 217-229.
             !
             ! A. Friend       : trans = 0.251 + 0.509 * (1.0 - cloud(i))
             ! Spitters et al. : trans = 0.200 + 0.560 * (1.0 - cloud(i))
             !
             ! we are using the values from A. Friend
             !
             trans = 0.251_r8 + 0.509_r8 * (1.0_r8 - cloud(i,jb)) 
             !
             ! calculate the fraction of indirect (diffuse) solar radiation
             ! based upon the cloud cover
             !
             ! note that these relationships typically are measured for either
             ! monthly or daily timescales, and may not be exactly appropriate
             ! for hourly calculations -- however, in ibis, cloud cover is fixed
             ! through the entire day so it may not make much difference
             !
             ! method i --
             !
             ! we use a simple empirical relationships from Nikolov and Zeller (1992)
             !
             ! Nikolov, N. and K.F. Zeller, 1992:  A solar radiation algorithm for ecosystem
             ! dynamics models, Ecological Modelling, 61, 149-168.
             !
             fdiffuse = 1.0045_r8 + 0.0435_r8 * trans - 3.5227_r8 * trans**2 + 2.6313_r8 * trans**3
             !
             IF (trans.GT.0.75_r8) fdiffuse = 0.166_r8
             !
             ! method ii --
             !
             ! another method was suggested by Spitters et al. (1986) based on
             ! long-term data from the Netherlands
             !
             ! Spitters et al., 1986: Separating the diffuse and direct component
             ! of global radiation and its implications for modeling canopy
             ! photosynthesis, Part I: Components of incoming radiation,
             ! Agricultural and Forest Meteorology, 38, 217-229.
             !
             !       if ((trans.eq.0.00).and.(trans.lt.0.07)) then
             !         fdiffuse = 1.0
             !       else if ((trans.ge.0.07).and.(trans.lt.0.35)) then
             !         fdiffuse = 1.0 - 2.3 * (trans - 0.07)**2
             !       else if ((trans.ge.0.35).and.(trans.lt.0.75)) then
             !         fdiffuse = 1.33 - 1.46 * trans
             !       else
             !         fdiffuse = 0.23
             !       endif
             !
             ! do for each waveband
             !

             !
             ! calculate the fraction in each waveband
             !
             !        frac = 0.46 + 0.08 * float(ib - 1)

             frac = 0.427_r8 + 0.146_r8 * REAL((ib - 1),KIND=r8)
             !
             ! calculate the direct and indirect solar radiation
             !
             solad(i,ib) = swdown(i,jb) * frac *  (1.0_r8 - fdiffuse)               !
             solai(i,ib) = swdown(i,jb) * frac *  fdiffuse

             cosz_sib   (i,jb) =coszen (i,jb)
             radvbc_sib (i,jb) =solad  (i,1) 
             radnbc_sib (i,jb) =solad  (i,2) 
             radvdc_sib (i,jb) =solai  (i,1) 
             radndc_sib (i,jb) =solai  (i,2) 
             !
          END DO
       END DO
    END DO

       DO jb=1,jbMax
          DO i=1,nCols
            swradtot(i,jb)    =radvbc_sib (i,jb)+ radnbc_sib (i,jb)+ radvdc_sib (i,jb)+ radndc_sib (i,jb)
            diffswradtot(i,jb)=swdown(i,jb)-swradtot(i,jb)
          END DO
       END DO

       DO jb=1,jbMax
          DO i=1,nCols
             radvbc_sib (i,jb) =radvbc_sib (i,jb)  + 0.25_r8*diffswradtot(i,jb)
             radnbc_sib (i,jb) =radnbc_sib (i,jb)  + 0.25_r8*diffswradtot(i,jb)
             radvdc_sib (i,jb) =radvdc_sib (i,jb)  + 0.25_r8*diffswradtot(i,jb)
             radndc_sib (i,jb) =radndc_sib (i,jb)  + 0.25_r8*diffswradtot(i,jb)
          END DO
       END DO

    RETURN
100 STOP 
  END SUBROUTINE RunAtmosModel1D


  SUBROUTINE RunAtmosModel2D(rec,jdt,tod,nband,dt,idate,idatec,idatep,DELTAIN,start,iMask)
    IMPLICIT NONE
    INTEGER       , INTENT(INOUT) :: rec
    REAL(KIND=r8), INTENT(IN   ) :: tod
    INTEGER       , INTENT(IN   ) :: jdt
    INTEGER       , INTENT(IN   ) :: nband  ! INTENT(IN    ) :: nband
    REAL(KIND=r8), INTENT(IN   ) :: dt     ! INTENT(IN    ) :: dt 
    INTEGER       , INTENT(IN   ) :: idate (4) ! INTENT(IN    ) :: idatec(:)
    INTEGER       , INTENT(IN   ) :: idatec(4) ! INTENT(IN    ) :: idatec(:)
    INTEGER       , INTENT(IN   ) :: idatep(4)
    REAL(KIND=r8), INTENT(IN   ) :: DELTAIN
    LOGICAL       , INTENT(IN   ) :: start
    INTEGER(KIND=i8)   , INTENT(IN   ) :: imask      (ibMax,  jbMax)
    INTEGER :: nCols
    INTEGER        :: jb
    INTEGER        :: ib
    INTEGER        :: i
    INTEGER        :: ios
    REAL(KIND=r8) :: rtime
    REAL(KIND=r8) :: time
    REAL(KIND=r8) :: time2
    REAL(KIND=r8) :: jday
    REAL(KIND=r8) :: orbit
    REAL(KIND=r8) :: angle
    REAL(KIND=r8) :: xdecl
    REAL(KIND=r8) :: frac
    REAL(KIND=r8) :: sw
    REAL(KIND=r8) :: year
    REAL(KIND=r8) :: decmax
    REAL(KIND=r8) :: sols
    REAL(KIND=r8) :: season
    REAL(KIND=r8) :: day
    REAL(KIND=r8) :: dec
    REAL(KIND=r8) :: rfd
    REAL(KIND=r8) :: sind
    REAL(KIND=r8) :: cosd
    REAL(KIND=r8) :: hac
    REAL(KIND=r8) :: h
    REAL(KIND=r8) :: sr
    REAL(KIND=r8) :: ss
    REAL(KIND=r8) :: dawn 
    REAL(KIND=r8) :: dusk  
    REAL(KIND=r8) :: coshr 
    REAL(KIND=r8) :: sunang
    REAL(KIND=r8) :: trans
    REAL(KIND=r8) :: btime 
    REAL(KIND=r8) :: ctime
    REAL(KIND=r8) :: atime 
    REAL(KIND=r8) :: fdiffuse
    REAL(KIND=r8) :: xlati (jbMax)
    REAL(KIND=r8) :: coszen(ibMax,jbMax)
    REAL(KIND=r8) :: cloud(ibMax,jbMax)
    REAL(KIND=r8) :: solad(ibMax,nband)
    REAL(KIND=r8) :: solai(ibMax,nband)
    REAL(KIND=r8) :: var(ibMax,jbMax)
    CHARACTER(LEN=255) :: namec
    CHARACTER(LEN=255) :: namep
    INTEGER            :: open_status
    INTEGER            :: reclen
    INTEGER            :: irec
    INTEGER            :: j
    INTEGER, PARAMETER   :: input_unit=200
    REAL(KIND=r8) , PARAMETER    :: yrl=365.2500_r8
    INTEGER, PARAMETER   :: ok = 0
    REAL(KIND=r8)   , PARAMETER :: alon   =  0.0_r8
    REAL (KIND=r8), PARAMETER   :: pie  =              3.1415926e0_r8! constant pi=3.1415926e0   
    REAL (KIND=r8), PARAMETER   :: pai12  =              pie/12.0_r8 ! pi/12

    LOGICAL :: here
    CHARACTER(LEN=255) :: PREFIX='GCM2IBIS'
    CHARACTER(LEN=255) :: SUFFIX='.dat'
    INTEGER           :: iyeard 
    INTEGER           :: imonth
    INTEGER           :: juliand
    REAL(KIND=r8)    :: frachours
    REAL(KIND=r8)    :: hours
    REAL(KIND=r8)    :: xua
    REAL(KIND=r8)    :: xta
    REAL(KIND=r8)    :: xprec
    REAL(KIND=r8)    :: xsin
    REAL(KIND=r8)    :: xlin
    REAL(KIND=r8)    :: xur
    REAL(KIND=r8)    :: xes
    REAL(KIND=r8)    :: xe
    REAL(KIND=r8)    :: xqa
    REAL(KIND=r8)    :: xpslc
    REAL(KIND=r8)    :: xlh=0
    REAL(KIND=r8)    :: xsh=0
    REAL(KIND=r8)    :: sdelt
    REAL(KIND=r8)    :: ratio
    REAL(KIND=r8)    :: etime
    REAL(KIND=r8)    :: xday
    REAL(KIND=r8)    :: sindel
    REAL(KIND=r8)    :: cosdel
    REAL(KIND=r8)    :: fimxi 
    REAL(KIND=r8)    :: frh
    REAL(KIND=r8)    :: cos2  (ibMax)
    REAL(KIND=r8)    :: zenith1 (ibMax)
    REAL(KIND=r8)    :: zenith2 (ibMax)
    REAL(KIND=r8)    :: lonrad(ibMax)

    REAL(KIND=r8)    :: xua_stepc(ibMax,jbMax)
    REAL(KIND=r8)    :: xta_stepc(ibMax,jbMax)
    REAL(KIND=r8)    :: xpres_stepc(ibMax,jbMax)
    REAL(KIND=r8)    :: xprec_stepc(ibMax,jbMax)
    REAL(KIND=r8)    :: xsin_stepc(ibMax,jbMax)
    REAL(KIND=r8)    :: xlin_stepc(ibMax,jbMax)
    REAL(KIND=r8)    :: xur_stepc(ibMax,jbMax)
    REAL(KIND=r8)    :: xqa_stepc(ibMax,jbMax)

    REAL(KIND=r8)    :: xua_stepp(ibMax,jbMax)
    REAL(KIND=r8)    :: xta_stepp(ibMax,jbMax)
    REAL(KIND=r8)    :: xpres_stepp(ibMax,jbMax)
    REAL(KIND=r8)    :: xprec_stepp(ibMax,jbMax)
    REAL(KIND=r8)    :: xsin_stepp(ibMax,jbMax)
    REAL(KIND=r8)    :: xlin_stepp(ibMax,jbMax)
    REAL(KIND=r8)    :: xur_stepp(ibMax,jbMax)
    REAL(KIND=r8)    :: xqa_stepp(ibMax,jbMax)
    REAL(KIND=r8)    :: swradtot    (ibMax,jbMax) 
    REAL(KIND=r8)    :: diffswradtot(ibMax,jbMax)
    
    open_status=ok
    nCols=ibMax
    time=tod
    IF(MOD(tod,DELTAIN) == 0.0_r8 .OR. jdt == 1)THEN
          step=0
          rec=rec+1
          IF(TRIM(icn_data) == 'ERA5')THEN
             !PK CALL RunReadNETCDF(rec,ncid_prec,latic,latip,longc,longp,xprec_stepc,xprec_stepp,nFile=1)
             !PK CALL RunReadNETCDF(rec,ncid_temp,latic,latip,longc,longp,xta_stepc  ,xta_stepp  ,nFile=2)
             !PK CALL RunReadNETCDF(rec,ncid_rhum,latic,latip,longc,longp,xur_stepc  ,xur_stepp  ,nFile=3)
             !PK CALL RunReadNETCDF(rec,ncid_wind,latic,latip,longc,longp,xua_stepc  ,xua_stepp  ,nFile=4)
             !PK CALL RunReadNETCDF(rec,ncid_swrd,latic,latip,longc,longp,xsin_stepc ,xsin_stepp ,nFile=5)
             !PK CALL RunReadNETCDF(rec,ncid_lwrd,latic,latip,longc,longp,xlin_stepc ,xlin_stepp, nFile=6)
          ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
             !PK CALL RunReadNETCDF(rec,ncid_pres,latic,latip,longc,longp,xpres_stepc,xpres_stepp,nFile=1)
             !PK CALL RunReadNETCDF(rec,ncid_temp,latic,latip,longc,longp,xta_stepc  ,xta_stepp  ,nFile=2)
             !PK CALL RunReadNETCDF(rec,ncid_umes,latic,latip,longc,longp,xqa_stepc  ,xqa_stepp  ,nFile=3)
             !PK CALL RunReadNETCDF(rec,ncid_wind,latic,latip,longc,longp,xua_stepc  ,xua_stepp  ,nFile=4)
             !PK CALL RunReadNETCDF(rec,ncid_prec,latic,latip,longc,longp,xprec_stepc,xprec_stepp,nFile=5)
             !PK CALL RunReadNETCDF(rec,ncid_swrd,latic,latip,longc,longp,xsin_stepc ,xsin_stepp ,nFile=6)
             !PK CALL RunReadNETCDF(rec,ncid_lwrd,latic,latip,longc,longp,xlin_stepc ,xlin_stepp, nFile=7)
          ELSE
             PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
             STOP       
          END IF


       DO j=1,jbMax
          DO i=1,ibMax
             !IF (iMask(i,jbMax+1-j) >= 1) THEN

           t02mc(i,j) = xta_stepc(i,j)+273.16_r8                     !K
           !xpslc      = 997.0_r8                                     !hPa
           IF(TRIM(icn_data) == 'ERA5')THEN
              xpslc      = (1.15_r8*287.16*(xta_stepc(i,j)+273.16_r8) )/100.0_r8    ! (rho*R*T)  hPa
              xes=6.112*exp(17.67*xta_stepc(i,j)/(243.5+xta_stepc(i,j)))!hPa
              xe=(xur_stepc(i,j)/100.0_r8)*xes                          !hPa
              xqa  =xe*622.0/(xpslc-xe)                                 !g/kg
              q02mc(i,j) = xqa/1000.0_r8                                !kg/kg
           ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
              xpslc      = (1.15_r8*287.16*(xta_stepc(i,j)+273.16_r8) )/100.0_r8    ! (rho*R*T)  hPa
              q02mc(i,j) =   xqa_stepc(i,j)
           ELSE
             PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
             STOP       
           END IF
           uvesc(i,j) = xua_stepc(i,j)                   !m/s
           vvesc(i,j) = xua_stepc(i,j)                                       !m/s
           IF(TRIM(icn_data) == 'ERA5')THEN
              pslcc(i,j) = xpslc                                        !hPa
           ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
	      IF(xpres_stepc(i,j) < 100)  THEN
                 pslcc(i,j)       =xpslc   
	      ELSE
                 pslcc(i,j)       = xpres_stepc(i,j)   
	      END IF
           ELSE
             PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
             STOP       
           END IF

           IF(TRIM(icn_data) == 'ERA5')THEN
              precc(i,j) = xprec_stepc(i,j)*86400/DELTAIN               !mm/dia
           ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
              precc(i,j) = xprec_stepc(i,j)  !*86400/DELTAIN               !mm/dia
           ELSE
             PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
             STOP       
           END IF
 
           ocisc(i,j) = xsin_stepc(i,j)                              !w/m2
           olisc(i,j) = xlin_stepc(i,j)                              !w/m2
         !  END IF
          END DO
       END DO 
       IF(.FALSE.)THEN
          WRITE(0,*)'latic',MAXVAL(latic),MINVAL(latic)
          WRITE(0,*)'longc',MAXVAL(longc),MINVAL(longc)
          WRITE(0,*)'t02mc',MAXVAL(t02mc),MINVAL(t02mc)
          WRITE(0,*)'q02mc',MAXVAL(q02mc),MINVAL(q02mc)
          WRITE(0,*)'uvesc',MAXVAL(uvesc),MINVAL(uvesc)
          WRITE(0,*)'vvesc',MAXVAL(vvesc),MINVAL(vvesc)
          WRITE(0,*)'pslcc',MAXVAL(pslcc),MINVAL(pslcc)
          WRITE(0,*)'precc',MAXVAL(precc),MINVAL(precc)
          WRITE(0,*)'ocisc',MAXVAL(ocisc),MINVAL(ocisc)
          WRITE(0,*)'olisc',MAXVAL(olisc),MINVAL(olisc)
          WRITE(0,*)'slhfc',MAXVAL(slhfc),MINVAL(slhfc)
          WRITE(0,*)'sshfc',MAXVAL(sshfc),MINVAL(sshfc)
       END IF
       
        DO j=1,jbMax
           DO i=1,ibMax
           !IF (iMask(i,jbMax+1-j) >= 1) THEN

              t02mp(i,j) = xta_stepp(i,j)+273.16_r8                              !K
              !xpslc      = 1000.0_r8                                             !hPa
              IF(TRIM(icn_data) == 'ERA5')THEN
                 xpslc      = (1.15_r8*287.16*(xta_stepp(i,j)+273.16_r8) )/100.0_r8    ! (rho*R*T)  hPa
                 xes        = 6.112*exp(17.67*xta_stepp(i,j)/(243.5+xta_stepp(i,j)))!hPa
                 xe         = (xur_stepp(i,j)/100.0_r8)*xes                         !hPa
                 xqa        = xe*622.0/(xpslc-xe)                                   !g/kg
                 q02mp(i,j) = xqa/1000.0_r8                                         !kg/kg
              ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
                 xpslc      = (1.15_r8*287.16*(xta_stepp(i,j)+273.16_r8) )/100.0_r8    ! (rho*R*T)  hPa
                 q02mp(i,j) =   xqa_stepp(i,j)
              ELSE
                 PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                 STOP       
              END IF
              uvesp(i,j) = xua_stepp(i,j)                            !m/s
              vvesp(i,j) = xua_stepp(i,j)                                               !m/s
              IF(TRIM(icn_data) == 'ERA5')THEN
                 pslcp(i,j) = xpslc                                        !hPa
              ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
	         IF(xpres_stepp(i,j) < 100)  THEN
                    pslcp(i,j)       =xpslc   
	         ELSE
                    pslcp(i,j)       = xpres_stepp(i,j)     
	         END IF

              ELSE
                PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                STOP       
              END IF
              IF(TRIM(icn_data) == 'ERA5')THEN
                 precp(i,j) = xprec_stepp(i,j)*86400/DELTAIN                        !mm/dia
              ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
                 precp(i,j) = xprec_stepp(i,j)!*86400/DELTAIN                        !mm/dia
              ELSE
                PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
                STOP       
              END IF

              ocisp(i,j) = xsin_stepp(i,j)                                       !w/m2
              olisp(i,j) = xlin_stepp(i,j)                                       !w/m2
              slhfp(i,j) = xlh
              sshfp(i,j) = xsh
          !   END IF
           END DO
        END DO     
        IF(.FALSE.)THEN
           WRITE(0,*)'latip',MAXVAL(latip),MINVAL(latip)
           WRITE(0,*)'longp',MAXVAL(longp),MINVAL(longp)
           WRITE(0,*)'t02mp',MAXVAL(t02mp),MINVAL(t02mp)
           WRITE(0,*)'q02mp',MAXVAL(q02mp),MINVAL(q02mp)
           WRITE(0,*)'uvesp',MAXVAL(uvesp),MINVAL(uvesp)
           WRITE(0,*)'vvesp',MAXVAL(vvesp),MINVAL(vvesp)
           WRITE(0,*)'pslcp',MAXVAL(pslcp),MINVAL(pslcp)
           WRITE(0,*)'precp',MAXVAL(precp),MINVAL(precp)
           WRITE(0,*)'ocisp',MAXVAL(ocisp),MINVAL(ocisp)
           WRITE(0,*)'olisp',MAXVAL(olisp),MINVAL(olisp)
           WRITE(0,*)'slhfp',MAXVAL(slhfp),MINVAL(slhfp)
           WRITE(0,*)'sshfp',MAXVAL(sshfp),MINVAL(sshfp)
        END IF
    END IF

    DO jb=1,jbMax
       DO ib=1,nCols
      !     IF (iMask(ib,jb) >= 1) THEN

          gu        (ib,jbMax+1-jb) =  0.00000_r8
          gv        (ib,jbMax+1-jb) =  0.00000_r8
          gt        (ib,jbMax+1-jb) =  0.00000_r8
          gq        (ib,jbMax+1-jb) =  0.00000_r8
          gps       (ib,jbMax+1-jb) =  0.00000_r8
          var       (ib,        jb) =  0.00000_r8
          dcupr_sib (ib,jbMax+1-jb) =  0.00000_r8
          dlspr_sib (ib,jbMax+1-jb) =  0.00000_r8
          swdown    (ib,jbMax+1-jb) =  0.00000_r8
          dlwbot_sib(ib,jbMax+1-jb) =  0.00000_r8
          sigki                     =  0.00000_r8
          delsig                    =  0.00000_r8
          tseam     (ib,jbMax+1-jb) =  0.00000_r8
          tsea      (ib,jbMax+1-jb) =  0.00000_r8
          slrad     (ib,jbMax+1-jb) =  0.00000_r8
          slhf      (ib,jbMax+1-jb) =  0.00000_r8
          sshf      (ib,jbMax+1-jb) =  0.00000_r8
          zorl      (ib,jbMax+1-jb) =  0.00000_r8
          tsurf     (ib,jbMax+1-jb) =  0.00000_r8
          qsurf     (ib,jbMax+1-jb) =  0.00000_r8
          tmtx      (ib,:,jbMax+1-jb) = 0.00000_r8
          qmtx      (ib,:,jbMax+1-jb) = 0.00000_r8
          umtx      (ib,:,jbMax+1-jb) = 0.00000_r8
       !  END IF    
       END DO     


        !pi --- 180
        ! y      x
        !
        ! y =  X*pi/180 !radianos
        ! 
        lati   (jbMax+1-jb)    = (latic(1,jb)+90.0_r8)*pi/180.0_r8
        xlati  (jb)            = (latic(1,jb)+90.0_r8)*pi/180.0_r8
       DO ib=1,nCols
          lonrad(ib)                = longc(ib,1)
           IF (iMask(ib,jbMax+1-jb) >= 1) THEN

          gu        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*uvesc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*uvesp(ib,jb)
          gv        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*vvesc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*vvesp(ib,jb)
          gt        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          gq        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*q02mp(ib,jb)
          gps       (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*pslcc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*pslcp(ib,jb)
          var       (ib,        jb) = ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*precp(ib,jb)
          IF(TRIM(icn_data) == 'ERA5')THEN
              dcupr_sib (ib,jbMax+1-jb) = (var(ib,jb)/86400.0_r8)*1.0e+3_r8
             !          PRINT*,  var       (ib,        jb), ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb), ((step*dt)/(DELTAIN))*precp(ib,jb)
              dlspr_sib (ib,jbMax+1-jb) = 0.0_r8
          ELSE IF(TRIM(icn_data) == 'GLDAS')THEN
              dcupr_sib (ib,jbMax+1-jb) = var(ib,jb)*1.0e+3_r8!/86400.0_r8)*1.0e+3_r8
!             PRINT*,  var       (ib,        jb), ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb), ((step*dt)/(DELTAIN))*precp(ib,jb)
              dlspr_sib (ib,jbMax+1-jb) = 0.0_r8
          ELSE
             PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS'
             STOP       
           END IF
          swdown    (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*ocisc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*ocisp(ib,jb)
          dlwbot_sib(ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*olisp(ib,jb)
          sigki                     = 1.0_r8
          delsig                    = 0.010000_r8
          tseam     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          tsea      (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          slrad     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*olisp(ib,jb)
          slhf      (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*slhfc(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*slhfp(ib,jb)
          sshf      (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*sshf(ib,jb) &
                                            + ((step*dt)/(DELTAIN))*sshfp(ib,jb)
    ! zorl........zorl (i)= 100.0 *zgrav*speedm(i)*rhi(i)
    !             zgrav =0.032 /grav

          zorl      (ib,jbMax+1-jb) = 100.0_r8*(0.32_r8/9.8_r8)*(sqrt( gu(ib,jbMax+1-jb)**2 + gv(ib,jbMax+1-jb)**2))*0.01
          tsurf     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
                                      + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          qsurf     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) &
                                      + ((step*dt)/(DELTAIN))*q02mp(ib,jb)
          tmtx      (ib,:,jbMax+1-jb) = 0.010000_r8
          qmtx      (ib,:,jbMax+1-jb) = 0.010000_r8
          umtx      (ib,:,jbMax+1-jb) = 0.010000_r8
         endif
       END DO
    END DO

    step=step+1
    
                IF(.TRUE.)THEN
                   WRITE(0,*)'latip',MAXVAL(latip),MINVAL(latip)
                   WRITE(0,*)'longp',MAXVAL(longp),MINVAL(longp)
                   WRITE(0,*)'t02mp',MAXVAL(gt),MINVAL(gt)
                   WRITE(0,*)'q02mp',MAXVAL(gq),MINVAL(gq)
                   WRITE(0,*)'uvesp',MAXVAL(gu),MINVAL(gu)
                   WRITE(0,*)'vvesp',MAXVAL(gv),MINVAL(gv)
                   WRITE(0,*)'pslcp',MAXVAL(gps),MINVAL(gps)
                   WRITE(0,*)'precp',MAXVAL(dcupr_sib),MINVAL(dcupr_sib)
                   WRITE(0,*)'ocisp',MAXVAL(swdown),MINVAL(swdown)
                   WRITE(0,*)'olisp',MAXVAL(dlwbot_sib),MINVAL(dlwbot_sib)
                   WRITE(0,*)'slhfp',MAXVAL(slhfp),MINVAL(slhfp)
                   WRITE(0,*)'sshfp',MAXVAL(sshfp),MINVAL(sshfp)
                END IF

    !
    !     mon is the month used for vegetation data input
    !
    DO jb=1,jbMax
       !DO ib=1,ibMax,1
       !   month(ib,jb)=idatec(2)
       !   IF((((lati(jb)*180.0_r8)/3.1415926e0_r8)-90.0_r8)  > 0.0_r8 ) THEN
       !      month(ib,jb)  =  month(ib,jb) + 6
       !      IF(month(ib,jb).GE.13) month(ib,jb) = month(ib,jb)-12
       !   END IF
       !END DO
       !
       !     computation of astronomical parameters
       !     sdelt ;solar inclination
       !     etime ;correction factor to local time
       !     ratio ;factor relating to the distance between the earth and the sun
       !
       CALL radtim(idatec,sdelt ,ratio ,etime ,tod   ,xday  ,yrl)
       sw=1370.0_r8*ratio

       sindel = SIN(sdelt)
       cosdel = COS(sdelt)
       fimxi  = 24.0e0_r8 /360.0_r8
       ctime  = alon/15.0e0_r8
       cos2   = 0.0e0_r8
       frh    = ( MOD(tod+0.03125_r8,3600.0_r8)-0.03125_r8)/3600.0_r8
       DO ib=1,ibMax
          zenith1(ib)  = sindel*COS(xlati(jb))
       ENDDO
       DO ib=1,ibMax
          btime       = fimxi*lonrad(ib)+ctime
          atime       = etime+pai12*(12.0_r8-idatec(1)-frh-btime)
          zenith2 (ib) = cosdel*SIN(xlati(jb))*COS(atime)
          coszen  (ib,jb) =  MAX (0.0_r8,zenith1(ib) + zenith2(ib))
          
          sunang = MAX( 0.01_r8, coszen  (ib,jb) )
          
          cloud(ib,jb) = (sw * sunang - swdown(ib,jb)) / (sw * sunang)
          cloud(ib,jb) = MAX(cloud(ib,jb),0.0_r8)
          cloud(ib,jb) = MIN(cloud(ib,jb),1.0_r8)
          !cloud(ib,jb) = MAX(0.58_r8     ,cloud(ib,jb))
       END DO

    END DO
    !
    ! do for each waveband
    !
    DO ib = 1, nband
       DO jb=1,jbMax
          DO i=1,nCols
             ! calculate the solar transmission through the atmosphere
             ! using simple linear function of tranmission and cloud cover
             !
             ! note that the 'cloud cover' data is typically obtained from
             ! sunshine hours -- not direct cloud observations
             !
             ! where, cloud cover = 1 - sunshine fraction 
             !
             ! different authors present different values for the slope and 
             ! intercept terms of this equation
             !
             ! Friend, A: Parameterization of a global daily weather generator for
             ! terrestrial ecosystem and biogeochemical modelling, Ecological 
             ! Modelling
             !
             ! Spitters et al., 1986: Separating the diffuse and direct component
             ! of global radiation and its implications for modeling canopy
             ! photosynthesis, Part I: Components of incoming radiation,
             ! Agricultural and Forest Meteorology, 38, 217-229.
             !
             ! A. Friend       : trans = 0.251 + 0.509 * (1.0 - cloud(i))
             ! Spitters et al. : trans = 0.200 + 0.560 * (1.0 - cloud(i))
             !
             ! we are using the values from A. Friend
             !
             trans = 0.251_r8 + 0.509_r8 * (1.0_r8 - cloud(i,jb)) 
             !
             ! calculate the fraction of indirect (diffuse) solar radiation
             ! based upon the cloud cover
             !
             ! note that these relationships typically are measured for either
             ! monthly or daily timescales, and may not be exactly appropriate
             ! for hourly calculations -- however, in ibis, cloud cover is fixed
             ! through the entire day so it may not make much difference
             !
             ! method i --
             !
             ! we use a simple empirical relationships from Nikolov and Zeller (1992)
             !
             ! Nikolov, N. and K.F. Zeller, 1992:  A solar radiation algorithm for ecosystem
             ! dynamics models, Ecological Modelling, 61, 149-168.
             !
             fdiffuse = 1.0045_r8 + 0.0435_r8 * trans - 3.5227_r8 * trans**2 + 2.6313_r8 * trans**3
             !
             IF (trans.GT.0.75_r8) fdiffuse = 0.166_r8
             !
             ! method ii --
             !
             ! another method was suggested by Spitters et al. (1986) based on
             ! long-term data from the Netherlands
             !
             ! Spitters et al., 1986: Separating the diffuse and direct component
             ! of global radiation and its implications for modeling canopy
             ! photosynthesis, Part I: Components of incoming radiation,
             ! Agricultural and Forest Meteorology, 38, 217-229.
             !
             !       if ((trans.eq.0.00).and.(trans.lt.0.07)) then
             !         fdiffuse = 1.0
             !       else if ((trans.ge.0.07).and.(trans.lt.0.35)) then
             !         fdiffuse = 1.0 - 2.3 * (trans - 0.07)**2
             !       else if ((trans.ge.0.35).and.(trans.lt.0.75)) then
             !         fdiffuse = 1.33 - 1.46 * trans
             !       else
             !         fdiffuse = 0.23
             !       endif
             !
             ! do for each waveband
             !

             !
             ! calculate the fraction in each waveband
             !
             !        frac = 0.46 + 0.08 * float(ib - 1)

             frac = 0.427_r8 + 0.146_r8 * REAL((ib - 1),KIND=r8)
             !
             ! calculate the direct and indirect solar radiation
             !
             solad(i,ib) = swdown(i,jb) * frac *  (1.0_r8 - fdiffuse)               !
             solai(i,ib) = swdown(i,jb) * frac *  fdiffuse

             cosz_sib   (i,jb) =coszen (i,jb)
             radvbc_sib (i,jb) =solad  (i,1) 
             radnbc_sib (i,jb) =solad  (i,2) 
             radvdc_sib (i,jb) =solai  (i,1) 
             radndc_sib (i,jb) =solai  (i,2) 
             !
          END DO
       END DO
    END DO

       DO jb=1,jbMax
          DO i=1,nCols
            swradtot(i,jb)    =radvbc_sib (i,jb)+ radnbc_sib (i,jb)+ radvdc_sib (i,jb)+ radndc_sib (i,jb)
            diffswradtot(i,jb)=swdown(i,jb)-swradtot(i,jb)
          END DO
       END DO

       DO jb=1,jbMax
          DO i=1,nCols
             radvbc_sib (i,jb) =radvbc_sib (i,jb)  + 0.25_r8*diffswradtot(i,jb)
             radnbc_sib (i,jb) =radnbc_sib (i,jb)  + 0.25_r8*diffswradtot(i,jb)
             radvdc_sib (i,jb) =radvdc_sib (i,jb)  + 0.25_r8*diffswradtot(i,jb)
             radndc_sib (i,jb) =radndc_sib (i,jb)  + 0.25_r8*diffswradtot(i,jb)
          END DO
       END DO

    RETURN
100 STOP 
  END SUBROUTINE RunAtmosModel2D





  SUBROUTINE RunAtmosModel(jdt,tod,nband,dt,idate,idatec,idatep,DELTAIN)
    IMPLICIT NONE
    REAL(KIND=r8), INTENT(IN   ) :: tod
    INTEGER       , INTENT(IN   ) :: jdt
    INTEGER       , INTENT(IN   ) :: nband  ! INTENT(IN    ) :: nband
    REAL(KIND=r8), INTENT(IN   ) :: dt     ! INTENT(IN    ) :: dt 
    INTEGER       , INTENT(IN   ) :: idate (4) ! INTENT(IN    ) :: idatec(:)
    INTEGER       , INTENT(IN   ) :: idatec(4) ! INTENT(IN    ) :: idatec(:)
    INTEGER       , INTENT(IN   ) :: idatep(4)
    REAL(KIND=r8), INTENT(IN   ) :: DELTAIN
    INTEGER :: nCols
    INTEGER        :: jb
    INTEGER        :: ib
    INTEGER        :: i
    INTEGER        :: ios
    REAL(KIND=r8) :: rtime
    REAL(KIND=r8) :: time
    REAL(KIND=r8) :: time2
    REAL(KIND=r8) :: jday
    REAL(KIND=r8) :: orbit
    REAL(KIND=r8) :: angle
    REAL(KIND=r8) :: xdecl
    REAL(KIND=r8) :: frac
    REAL(KIND=r8) :: sw
    REAL(KIND=r8) :: year
    REAL(KIND=r8) :: decmax
    REAL(KIND=r8) :: sols
    REAL(KIND=r8) :: season
    REAL(KIND=r8) :: day
    REAL(KIND=r8) :: dec
    REAL(KIND=r8) :: rfd
    REAL(KIND=r8) :: sind
    REAL(KIND=r8) :: cosd
    REAL(KIND=r8) :: hac
    REAL(KIND=r8) :: h
    REAL(KIND=r8) :: sr
    REAL(KIND=r8) :: ss
    REAL(KIND=r8) :: dawn 
    REAL(KIND=r8) :: dusk  
    REAL(KIND=r8) :: coshr 
    REAL(KIND=r8) :: sunang
    REAL(KIND=r8) :: trans
    REAL(KIND=r8) :: btime 
    REAL(KIND=r8) :: ctime
    REAL(KIND=r8) :: atime 
    REAL(KIND=r8) :: fdiffuse
    REAL(KIND=r8) :: xlati (jbMax)
    REAL(KIND=r8) :: coszen(ibMax,jbMax)
    REAL(KIND=r8) :: cloud(ibMax,jbMax)
    REAL(KIND=r8) :: swdown(ibMax,jbMax)
    REAL(KIND=r8) :: solad(ibMax,nband)
    REAL(KIND=r8) :: solai(ibMax,nband)
    REAL(KIND=r8) :: var(ibMax,jbMax)
    CHARACTER(LEN=255) :: namec
    CHARACTER(LEN=255) :: namep
    INTEGER            :: open_status
    INTEGER            :: reclen
    INTEGER            :: irec
    INTEGER, PARAMETER   :: input_unit=200
    REAL(KIND=r8) , PARAMETER    :: yrl=365.2500_r8
    INTEGER, PARAMETER   :: ok = 0
    REAL(KIND=r8)   , PARAMETER :: alon   =  0.0_r8
    REAL (KIND=r8), PARAMETER   :: pie  =              3.1415926e0_r8! constant pi=3.1415926e0   
    REAL (KIND=r8), PARAMETER   :: pai12  =              pie/12.0_r8 ! pi/12

    LOGICAL :: here
    CHARACTER(LEN=255) :: PREFIX='GCM2IBIS'
    CHARACTER(LEN=255) :: SUFFIX='.dat'
    REAL(KIND=r8)    :: sdelt
    REAL(KIND=r8)    :: ratio
    REAL(KIND=r8)    :: etime
    REAL(KIND=r8)    :: xday
    REAL(KIND=r8)    :: sindel
    REAL(KIND=r8)    :: cosdel
    REAL(KIND=r8)    :: fimxi 
    REAL(KIND=r8)    :: frh
    REAL(KIND=r8)    :: cos2  (ibMax)
    REAL(KIND=r8)    :: zenith1 (ibMax)
    REAL(KIND=r8)    :: zenith2 (ibMax)
    REAL(KIND=r8)    :: lonrad(ibMax)
    nCols=ibMax
    time=tod
    IF(MOD(tod,DELTAIN) == 0.0_r8 .OR. jdt == 1)THEN
       step=0
       WRITE(namec,'(A8,2(I4.4,3i2.2),A4)')TRIM(PREFIX),&
       idate(4),idate(2),idate(3),idate(1),&
       idatec(4),idatec(2),idatec(3),idatec(1),&
       TRIM(SUFFIX)    
       INQUIRE (FILE='./inputdata/'//TRIM(namec),EXIST=here)
       IF (.NOT.here) THEN
          WRITE(0,*)'FILE ','./inputdata/'//TRIM(namec), 'Not Exist';STOP
       ELSE
          INQUIRE(IOLENGTH=reclen)latic
          OPEN(input_unit,FILE='./inputdata/'//TRIM(namec),ACCESS='DIRECT',&
             FORM='UNFORMATTED',STATUS='OLD',RECL=reclen, &
             IOSTAT=open_status)
          IF( open_status == ok ) THEN
             irec=1
             READ(input_unit,rec=irec,iostat = ios)latic;IF(ios /= 0) STOP
             WRITE(0,*)'latic',MAXVAL(latic),MINVAL(latic)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)longc;IF(ios /= 0) STOP
             WRITE(0,*)'longc',MAXVAL(longc),MINVAL(longc)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)t02mc;IF(ios /= 0) STOP
             WRITE(0,*)'t02mc',MAXVAL(t02mc),MINVAL(t02mc)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)q02mc;IF(ios /= 0) STOP
             WRITE(0,*)'q02mc',MAXVAL(q02mc),MINVAL(q02mc)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)uvesc;IF(ios /= 0) STOP
             WRITE(0,*)'uvesc',MAXVAL(uvesc),MINVAL(uvesc)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)vvesc;IF(ios /= 0) STOP
             WRITE(0,*)'vvesc',MAXVAL(vvesc),MINVAL(vvesc)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)pslcc;IF(ios /= 0) STOP
             WRITE(0,*)'pslcc',MAXVAL(pslcc),MINVAL(pslcc)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)precc;IF(ios /= 0) STOP
             WRITE(0,*)'precc',MAXVAL(precc),MINVAL(precc)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)ocisc;IF(ios /= 0) STOP
             WRITE(0,*)'ocisc',MAXVAL(ocisc),MINVAL(ocisc)
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)olisc;IF(ios /= 0) STOP
             WRITE(0,*)'olisc',MAXVAL(olisc),MINVAL(olisc)
             CLOSE(input_unit,STATUS='KEEP')             
          ELSE
             WRITE(0,*)"Unable to OPEN file";STOP
          END IF
       END IF

       WRITE(namep,'(A8,2(I4.4,3i2.2),A4)')TRIM(PREFIX),&
       idate(4),idate(2),idate(3),idate(1),&
       idatep(4),idatep(2),idatep(3),idatep(1),&
       TRIM(SUFFIX)    
       
       INQUIRE (FILE='./inputdata/'//TRIM(namep),EXIST=here)
       IF (.NOT.here) THEN
          WRITE(0,*)'FILE ','./inputdata/'//TRIM(namep), 'Not Exist';STOP
       ELSE
          INQUIRE(IOLENGTH=reclen)latic
          OPEN(input_unit,FILE='./inputdata/'//TRIM(namep),ACCESS='DIRECT',&
             FORM='UNFORMATTED',STATUS='OLD',RECL=reclen, &
             IOSTAT=open_status)
          IF( open_status == ok ) THEN
             irec=1
             READ(input_unit,rec=irec,iostat = ios)latip;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)longp;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)t02mp;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)q02mp;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)uvesp;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)vvesp;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)pslcp;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)precp;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)ocisp;IF(ios /= 0) STOP
             irec=irec+1
             READ(input_unit,rec=irec,iostat = ios)olisp;IF(ios /= 0) STOP
             CLOSE(input_unit,STATUS='KEEP')             
          ELSE
             WRITE(0,*)"Unable to OPEN file";STOP
          END IF
       END IF 
       WRITE(0,*)'READ ',TRIM(namec),' and ',TRIM(namep)
    END IF
    DO jb=1,jbMax
        !pi --- 180
        ! y      x
        !
        ! y =  X*pi/180 !radianos
        ! 
        lati   (jbMax+1-jb)    = (latic(1,jb)+90.0_r8)*pi/180.0_r8
        xlati  (jbMax+1-jb)    = (latic(1,jb))*pi/180.0_r8
       DO ib=1,nCols
          lonrad(ib)                = longc(ib,1)
          gu        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*uvesc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*uvesp(ib,jb)
          gv        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*vvesc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*vvesp(ib,jb)
          gt        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          gq        (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*q02mp(ib,jb)
          gps       (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*pslcc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*pslcp(ib,jb)
          var       (ib,        jb) = ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*precp(ib,jb)
          dcupr_sib (ib,jbMax+1-jb) = (var(ib,jb)/86400.0_r8)*1.0e+3_r8
          dlspr_sib (ib,jbMax+1-jb) = 0.0_r8
          swdown    (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*ocisc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*ocisp(ib,jb)
          dlwbot_sib(ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*olisp(ib,jb)
          sigki                     = 1.0_r8
          delsig                    = 0.010000_r8
          tseam     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          tsea      (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          slrad     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*olisp(ib,jb)
    ! zorl........zorl (i)= 100.0 *zgrav*speedm(i)*rhi(i)
    !             zgrav =0.032 /grav

          zorl      (ib,jbMax+1-jb) = 100.0_r8*(0.32_r8/9.8_r8)*(sqrt( gu(ib,jbMax+1-jb)**2 + gv(ib,jbMax+1-jb)**2))*0.01
          tsurf     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*t02mp(ib,jb)
          qsurf     (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) &
                                    + ((step*dt)/(DELTAIN))*q02mp(ib,jb)
          tmtx      (ib,:,jbMax+1-jb) = 0.010000_r8
          qmtx      (ib,:,jbMax+1-jb) = 0.010000_r8
          umtx      (ib,:,jbMax+1-jb) = 0.010000_r8
       END DO
    END DO

    step=step+1
    

    !
    !     mon is the month used for vegetation data input
    !
    DO jb=1,jbMax
       !DO ib=1,ibMax,1
       !   month(ib,jb)=idatec(2)
       !   IF((((lati(jb)*180.0_r8)/3.1415926e0_r8)-90.0_r8)  > 0.0_r8 ) THEN
       !      month(ib,jb)  =  month(ib,jb) + 6
       !      IF(month(ib,jb).GE.13) month(ib,jb) = month(ib,jb)-12
       !   END IF
       !END DO
       !
       !     computation of astronomical parameters
       !     sdelt ;solar inclination
       !     etime ;correction factor to local time
       !     ratio ;factor relating to the distance between the earth and the sun
       !
       CALL radtim(idatec,sdelt ,ratio ,etime ,tod   ,xday  ,yrl)
       sw=1370.0_r8*ratio

       sindel = SIN(sdelt)
       cosdel = COS(sdelt)
       fimxi  = 24.0e0_r8 /360.0_r8
       ctime  = alon/15.0e0_r8
       cos2   = 0.0e0_r8
       frh    = ( MOD(tod+0.03125_r8,3600.0_r8)-0.03125_r8)/3600.0_r8
       DO ib=1,ibMax
          zenith1(ib)  = sindel*COS(lati(jb))
       ENDDO
       DO ib=1,ibMax
          btime       = fimxi*lonrad(ib)+ctime
          atime       = etime+pai12*(12.0_r8-idatec(1)-frh-btime)
          zenith2 (ib) = cosdel*SIN(lati(jb))*COS(atime)
          coszen  (ib,jb) =  MAX (0.0_r8,zenith1(ib) + zenith2(ib))
          
          sunang = MAX( 0.01_r8, coszen  (ib,jb) )
          
          cloud(ib,jb) = (sw * sunang - swdown(ib,jb)) / (sw * sunang)
          cloud(ib,jb) = MAX(cloud(ib,jb),0.0_r8)
          cloud(ib,jb) = MIN(cloud(ib,jb),1.0_r8)
          !cloud(ib,jb) = MAX(0.58_r8     ,cloud(ib,jb))
       END DO

    END DO
    !
    ! do for each waveband
    !
    DO ib = 1, nband
       DO jb=1,jbMax
          DO i=1,nCols
             ! calculate the solar transmission through the atmosphere
             ! using simple linear function of tranmission and cloud cover
             !
             ! note that the 'cloud cover' data is typically obtained from
             ! sunshine hours -- not direct cloud observations
             !
             ! where, cloud cover = 1 - sunshine fraction 
             !
             ! different authors present different values for the slope and 
             ! intercept terms of this equation
             !
             ! Friend, A: Parameterization of a global daily weather generator for
             ! terrestrial ecosystem and biogeochemical modelling, Ecological 
             ! Modelling
             !
             ! Spitters et al., 1986: Separating the diffuse and direct component
             ! of global radiation and its implications for modeling canopy
             ! photosynthesis, Part I: Components of incoming radiation,
             ! Agricultural and Forest Meteorology, 38, 217-229.
             !
             ! A. Friend       : trans = 0.251 + 0.509 * (1.0 - cloud(i))
             ! Spitters et al. : trans = 0.200 + 0.560 * (1.0 - cloud(i))
             !
             ! we are using the values from A. Friend
             !
             trans = 0.251_r8 + 0.509_r8 * (1.0_r8 - cloud(i,jb)) 
             !
             ! calculate the fraction of indirect (diffuse) solar radiation
             ! based upon the cloud cover
             !
             ! note that these relationships typically are measured for either
             ! monthly or daily timescales, and may not be exactly appropriate
             ! for hourly calculations -- however, in ibis, cloud cover is fixed
             ! through the entire day so it may not make much difference
             !
             ! method i --
             !
             ! we use a simple empirical relationships from Nikolov and Zeller (1992)
             !
             ! Nikolov, N. and K.F. Zeller, 1992:  A solar radiation algorithm for ecosystem
             ! dynamics models, Ecological Modelling, 61, 149-168.
             !
             fdiffuse = 1.0045_r8 + 0.0435_r8 * trans - 3.5227_r8 * trans**2 + 2.6313_r8 * trans**3
             !
             IF (trans.GT.0.75_r8) fdiffuse = 0.166_r8
             !
             ! method ii --
             !
             ! another method was suggested by Spitters et al. (1986) based on
             ! long-term data from the Netherlands
             !
             ! Spitters et al., 1986: Separating the diffuse and direct component
             ! of global radiation and its implications for modeling canopy
             ! photosynthesis, Part I: Components of incoming radiation,
             ! Agricultural and Forest Meteorology, 38, 217-229.
             !
             !       if ((trans.eq.0.00).and.(trans.lt.0.07)) then
             !         fdiffuse = 1.0
             !       else if ((trans.ge.0.07).and.(trans.lt.0.35)) then
             !         fdiffuse = 1.0 - 2.3 * (trans - 0.07)**2
             !       else if ((trans.ge.0.35).and.(trans.lt.0.75)) then
             !         fdiffuse = 1.33 - 1.46 * trans
             !       else
             !         fdiffuse = 0.23
             !       endif
             !
             ! do for each waveband
             !

             !
             ! calculate the fraction in each waveband
             !
             !        frac = 0.46 + 0.08 * float(ib - 1)

             frac = 0.827_r8 + 0.146_r8 * REAL((ib - 1),KIND=r8)
             !
             ! calculate the direct and indirect solar radiation
             !
             solad(i,ib) = sw * coszen(i,jb) * frac * trans * (1.0_r8 - fdiffuse)               !
             solai(i,ib) = sw * coszen(i,jb) * frac * trans * fdiffuse

             cosz_sib   (i,jb) =coszen (i,jb)
             radvbc_sib (i,jb) =solad  (i,1) 
             radnbc_sib (i,jb) =solad  (i,2) 
             radvdc_sib (i,jb) =solai  (i,1) 
             radndc_sib (i,jb) =solai  (i,2) 
             !
          END DO
       END DO
    END DO


    RETURN
100 STOP 
  END SUBROUTINE RunAtmosModel
  ! radtim :calculates the astronomical parameters: solar inclination,
  !         correction factor to local time, factor relating to the distance
  !         between earth and sun, and the julian day.
  SUBROUTINE radtim (id    ,delta ,ratio ,etime ,tod   ,xday  ,yrl)
    !
    !==========================================================================
    !
    !==========================================================================
    !  id(1)....hour(00/12)
    !  id(2)....month
    !  id(3)....day of month
    !  id(4)....year
    !  delta....solar inclination
    !  ratio....factor relating to the distance between the earth and the sun
    !  etime....correction factor to local time
    !  tod......model forecast time of day in seconds
    !  xday.....is julian day - 1 with fraction of day
    !  pai......constant pi=3.1415926
    !  yrl......length of year in days
    !  monl.....length of each month in days
    !==========================================================================
    INTEGER, INTENT(in ) :: id(4)
    REAL(KIND=r8),    INTENT(out) :: delta
    REAL(KIND=r8),    INTENT(out) :: ratio
    REAL(KIND=r8),    INTENT(out) :: etime
    REAL(KIND=r8),    INTENT(in ) :: tod
    REAL(KIND=r8),    INTENT(out) :: xday
    REAL(KIND=r8),    INTENT(in ) :: yrl

    REAL(KIND=r8),    PARAMETER :: day0=-1.0_r8
    REAL(KIND=r8),    PARAMETER :: f3600=3.6e3_r8
    REAL(KIND=r8)          :: psi

    !
    !     id is now assumed to be the current date and hour
    !
    xday=id(1)*f3600
    xday=xday+MOD(tod,f3600)
    xday=monday(id(2))+id(3)+xday/86400.0_r8

    IF (yrl == 365.25e0_r8) THEN
       xday=xday-MOD(id(4)+3,4)*0.25_r8
       IF(MOD(id(4),4).EQ.0.AND.id(2).GT.2)xday=xday+1.0e0_r8
    END IF

    xday= MOD(xday-1.0e0_r8,yrl)

    IF (xday > day0) THEN
       psi=2.0e0_r8*pai*xday/yrl
       delta=0.006918e0_r8-0.399912e0_r8*COS(   psi)+0.070257e0_r8*SIN(psi) &
            -0.006758e0_r8*COS(2.0e0_r8*psi)+0.000907e0_r8*SIN(2.0e0_r8*psi) &
            -0.002697e0_r8*COS(3.0e0_r8*psi)+0.001480e0_r8*SIN(3.0e0_r8*psi)
       ratio=1.000110e0_r8+0.034221e0_r8*COS(   psi)+0.001280e0_r8*SIN(psi) &
            +0.000719e0_r8*COS(2.0e0_r8*psi)+0.000077e0_r8*SIN(2.0e0_r8*psi)
       etime=0.000075e0_r8+0.001868e0_r8*COS(   psi)-0.032077e0_r8*SIN(psi) &
            -0.014615e0_r8*COS(2.0e0_r8*psi)-0.040849e0_r8*SIN(2.0e0_r8*psi)
    ELSE
       WRITE(0,20)id,tod,xday
       WRITE(0,20)id,tod,xday
       STOP 2020
    END IF

20  FORMAT(' BAD DATE IN RADTIM.  ID=',4I5,' TOD=',G16.8,' XDAY=', &
         G16.8)
  END SUBROUTINE radtim

  SUBROUTINE InitRadtim(monl)
    INTEGER, INTENT(in ) :: monl(12)
    INTEGER              :: m
    monday(1)=0
    DO m=2,12
       monday(m)=monday(m-1)+monl(m-1)
    END DO
  END SUBROUTINE InitRadtim

END MODULE AtmosModel
