New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15422 for NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DOM/dtatsd.F90 – NEMO

Ignore:
Timestamp:
2021-10-21T11:19:25+02:00 (3 years ago)
Author:
jcastill
Message:

Changes tested so that they can merged with the CO9 Met Office branch - jpmax_harmo should be 34 with FES14 tides, but the last components are not used anyway

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DOM/dtatsd.F90

    r14075 r15422  
    2121   ! 
    2222   USE in_out_manager  ! I/O manager 
     23   USE iom             ! IOM library 
    2324   USE lib_mpp         ! MPP library 
    2425 
     
    3132   !                                  !!* namtsd  namelist : Temperature & Salinity Data * 
    3233   LOGICAL , PUBLIC ::   ln_tsd_init   !: T & S data flag 
     34   LOGICAL , PUBLIC ::   ln_tsd_interp !: vertical interpolation flag 
    3335   LOGICAL , PUBLIC ::   ln_tsd_dmp    !: internal damping toward input data flag 
    3436 
    3537   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read) 
     38   INTEGER                              ::   jpk_init, inum_dta   
     39   INTEGER                              ::   id, linum   ! local integers   
     40   INTEGER                              ::   zdim(4) 
    3641 
    3742   !!---------------------------------------------------------------------- 
     
    5358      LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used 
    5459      ! 
    55       INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3   ! local integers 
     60      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5   ! local integers 
    5661      !! 
    5762      CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files 
    58       TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read 
    59       TYPE(FLD_N)                   ::   sn_tem, sn_sal 
     63      TYPE(FLD_N), DIMENSION(jpts+2)::   slf_i           ! array of namelist informations on the fields to read   
     64      TYPE(FLD_N)                   ::   sn_tem, sn_sal, sn_dep, sn_msk  
    6065      !! 
    6166      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_dmp, cn_dir, sn_tem, sn_sal 
     
    6368      ! 
    6469      !  Initialisation 
    65       ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0 
     70      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  ; ierr4 = 0  ;  ierr5 = 0 
    6671      ! 
    6772      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :  
     
    8085         WRITE(numout,*) '~~~~~~~~~~~~ ' 
    8186         WRITE(numout,*) '   Namelist namtsd' 
    82          WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init = ', ln_tsd_init 
    83          WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_dmp  = ', ln_tsd_dmp 
     87         WRITE(numout,*) '      Initialisation of ocean T & S with T & S input data  ln_tsd_init   = ', ln_tsd_init  
     88         WRITE(numout,*) '      Interpolation of initial conditions in the vertical  ln_tsd_interp = ', ln_tsd_interp  
     89         WRITE(numout,*) '      damping of ocean T & S toward T & S input data       ln_tsd_dmp    = ', ln_tsd_dmp 
    8490         WRITE(numout,*) 
    8591         IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_dmp ) THEN 
     
    94100         ln_tsd_init = .FALSE. 
    95101      ENDIF 
     102      IF( ln_tsd_interp .AND. ln_tsd_dmp ) THEN   
     103            CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' )   ;   RETURN   
     104      ENDIF   
     105      IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN   
     106            CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' )   ;   RETURN   
     107      ENDIF 
    96108      ! 
    97109      !                             ! allocate the arrays (if necessary) 
    98110      IF( ln_tsd_init .OR. ln_tsd_dmp ) THEN 
    99111         ! 
    100          ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 
     112         IF( ln_tsd_interp ) THEN   
     113           ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information   
     114         ELSE   
     115           ALLOCATE( sf_tsd(jpts  ), STAT=ierr0 )    
     116         ENDIF 
    101117         IF( ierr0 > 0 ) THEN 
    102118            CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN 
    103119         ENDIF 
    104120         ! 
    105                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 ) 
    106          IF( sn_tem%ln_tint )   ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 
    107                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 ) 
    108          IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 
    109          ! 
    110          IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN 
     121         IF( ln_tsd_interp ) THEN   
     122            CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta )    
     123            id = iom_varid( inum_dta, sn_dep%clvar, zdim )   
     124            jpk_init = zdim(3)   
     125            IF(lwp) WRITE(numout,*) 'Dimension of vertical coordinate in ICs: ', jpk_init   
     126            CALL iom_close( inum_dta )   ! Close the input file   
     127            !   
     128                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr0 )   
     129            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 )   
     130                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr2 )   
     131            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 )     
     132                                 ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr4 )   
     133                                 ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr5 )   
     134         ELSE   
     135                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )   
     136            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )   
     137                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )   
     138            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )     
     139         ENDIF ! ln_tsd_interp   
     140         !   
     141         IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN  
    111142            CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN 
    112143         ENDIF 
    113144         !                         ! fill sf_tsd with sn_tem & sn_sal and control print 
    114145         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal 
     146         IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep   ;   slf_i(jp_msk) = sn_msk 
    115147         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print ) 
    116148         ! 
     
    138170      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data 
    139171      ! 
    140       INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies 
    141       INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers 
     172      INTEGER ::   ji, jj, jk, jl, jk_init   ! dummy loop indices   
     173      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1        ! local integers  
    142174      REAL(wp)::   zl, zi                             ! local scalars 
    143       REAL(wp), DIMENSION(jpk) ::  ztp, zsp   ! 1D workspace 
    144175      !!---------------------------------------------------------------------- 
    145176      ! 
     
    176207!!gm end 
    177208      ! 
    178       ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask 
    179       ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)  
    180       ! 
    181       IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==! 
    182          ! 
    183          IF( kt == nit000 .AND. lwp )THEN 
    184             WRITE(numout,*) 
    185             WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh' 
    186          ENDIF 
    187          ! 
    188          DO jj = 1, jpj                         ! vertical interpolation of T & S 
    189             DO ji = 1, jpi 
    190                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     209      IF( kt == nit000 .AND. lwp )THEN   
     210         WRITE(numout,*)   
     211         WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh'   
     212      ENDIF   
     213      !   
     214      IF( ln_tsd_interp ) THEN                 ! probably should use pointers in the following to make more readable   
     215      !   
     216         DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points   
     217            DO jj= 1, jpj   
     218               DO ji= 1, jpi 
    191219                  zl = gdept_0(ji,jj,jk) 
    192                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
    193                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
    194                      zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
    195                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
    196                      ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
    197                      zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
    198                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    199                      DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
    200                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    201                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    202                            ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
    203                            zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
     220                  IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN                     ! above the first level of data   
     221                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1)    
     222                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1)   
     223                  ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN          ! below the last level of data   
     224                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init)   
     225                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init)   
     226                  ELSE                                                             ! inbetween : vertical interpolation between jk_init & jk_init+1   
     227                     DO jk_init = 1, jpk_init-1                                    ! when  gdept(jk_init) < zl < gdept(jk_init+1)   
     228                        IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN       ! if there is no data fill down   
     229                           sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init)   
     230                           sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init)   
     231                        ENDIF   
     232                        IF((zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN   
     233                           zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / &   
     234                        &       (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init))   
     235                           ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + &   
     236                        &                         (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi   
     237                           ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + &   
     238                        &                         (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi 
    204239                        ENDIF 
    205240                     END DO 
    206241                  ENDIF 
    207242               END DO 
    208                DO jk = 1, jpkm1 
    209                   ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
    210                   ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
    211                END DO 
    212                ptsd(ji,jj,jpk,jp_tem) = 0._wp 
    213                ptsd(ji,jj,jpk,jp_sal) = 0._wp 
    214243            END DO 
    215244         END DO 
    216245         !  
     246         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:)   
     247         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:) 
    217248      ELSE                                !==   z- or zps- coordinate   ==! 
    218249         !                              
    219          ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask 
    220          ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:) 
     250         ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:) * tmask(:,:,:)  ! Mask   
     251         ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) * tmask(:,:,:) 
    221252         ! 
    222253         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
     
    248279                                        DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure 
    249280         IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta ) 
     281         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_dep)%fnow )     ! T arrays in the structure   
     282         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_msk)%fnow )     ! T arrays in the structure 
    250283                                        DEALLOCATE( sf_tsd              )     ! the structure itself 
    251284      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.