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 4147 for branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 – NEMO

Ignore:
Timestamp:
2013-11-04T12:51:55+01:00 (10 years ago)
Author:
cetlod
Message:

merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90

    r3294 r4147  
    160160      IF(lwp) WRITE(numout,*) '        orca ocean configuration' 
    161161 
    162 #if defined key_antarctic 
    163 #     include "ldfdyn_antarctic.h90" 
    164 #elif defined key_arctic 
    165 #     include "ldfdyn_arctic.h90" 
    166 #else 
    167       ! Read 2d integer array to specify western boundary increase in the 
    168       ! ===================== equatorial strip (20N-20S) defined at t-points 
    169  
    170       CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
    171       READ(inum,9101) clexp, iim, ijm 
    172       READ(inum,'(/)') 
    173       ifreq = 40 
    174       il1 = 1 
    175       DO jn = 1, jpidta/ifreq+1 
     162      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN 
     163! 
     164! 1.2 Modify ahm 
     165! -------------- 
     166         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean' 
     167         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm' 
     168         IF(lwp)WRITE(numout,*) '         north boundary increase' 
     169 
     170         ahm1(:,:) = ahm0 
     171         ahm2(:,:) = ahm0 
     172 
     173         ijpt0=max(1,min(49 -njmpp+1,jpj)) 
     174         ijpt1=max(0,min(49-njmpp+1,jpj-1)) 
     175         DO jj=ijpt0,ijpt1 
     176            ahm2(:,jj)=ahm0*2. 
     177            ahm1(:,jj)=ahm0*2. 
     178         END DO 
     179         ijpt0=max(1,min(48 -njmpp+1,jpj)) 
     180         ijpt1=max(0,min(48-njmpp+1,jpj-1)) 
     181         DO jj=ijpt0,ijpt1 
     182            ahm2(:,jj)=ahm0*1.9 
     183            ahm1(:,jj)=ahm0*1.75 
     184         END DO 
     185         ijpt0=max(1,min(47 -njmpp+1,jpj)) 
     186         ijpt1=max(0,min(47-njmpp+1,jpj-1)) 
     187         DO jj=ijpt0,ijpt1 
     188            ahm2(:,jj)=ahm0*1.5 
     189            ahm1(:,jj)=ahm0*1.25 
     190         END DO 
     191         ijpt0=max(1,min(46 -njmpp+1,jpj)) 
     192         ijpt1=max(0,min(46-njmpp+1,jpj-1)) 
     193         DO jj=ijpt0,ijpt1 
     194            ahm2(:,jj)=ahm0*1.1 
     195         END DO 
     196 
     197      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN 
     198! 1.2 Modify ahm  
     199! -------------- 
     200         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean' 
     201         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm' 
     202         IF(lwp)WRITE(numout,*) '         south and west boundary increase' 
     203 
     204 
     205         ahm1(:,:) = ahm0 
     206         ahm2(:,:) = ahm0 
     207 
     208         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj)) 
     209         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1)) 
     210         DO jj=ijpt0,ijpt1 
     211            ahm2(:,jj)=ahm0*2. 
     212            ahm1(:,jj)=ahm0*2. 
     213         END DO 
     214         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj)) 
     215         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1)) 
     216         DO jj=ijpt0,ijpt1 
     217            ahm2(:,jj)=ahm0*1.9 
     218            ahm1(:,jj)=ahm0*1.75 
     219         END DO 
     220         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj)) 
     221         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1)) 
     222         DO jj=ijpt0,ijpt1 
     223            ahm2(:,jj)=ahm0*1.5 
     224            ahm1(:,jj)=ahm0*1.25 
     225         END DO 
     226         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj)) 
     227         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1)) 
     228         DO jj=ijpt0,ijpt1 
     229            ahm2(:,jj)=ahm0*1.1 
     230         END DO 
     231      ELSE 
     232         ! Read 2d integer array to specify western boundary increase in the 
     233         ! ===================== equatorial strip (20N-20S) defined at t-points 
     234          
     235         CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
     236         READ(inum,9101) clexp, iim, ijm 
    176237         READ(inum,'(/)') 
    177          il2 = MIN( jpidta, il1+ifreq-1 ) 
    178          READ(inum,9201) ( ii, ji = il1, il2, 5 ) 
    179          READ(inum,'(/)') 
    180          DO jj = jpjdta, 1, -1 
    181             READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 ) 
    182          END DO 
    183          il1 = il1 + ifreq 
    184       END DO 
    185        
    186       DO jj = 1, nlcj 
    187          DO ji = 1, nlci 
    188             icof(ji,jj) = idata( mig(ji), mjg(jj) ) 
    189          END DO 
    190       END DO 
    191       DO jj = nlcj+1, jpj 
    192          DO ji = 1, nlci 
    193             icof(ji,jj) = icof(ji,nlcj) 
    194          END DO 
    195       END DO 
    196       DO jj = 1, jpj 
    197          DO ji = nlci+1, jpi 
    198             icof(ji,jj) = icof(nlci,jj) 
    199          END DO 
    200       END DO 
    201  
    202  9101 FORMAT(1x,a15,2i8) 
    203  9201 FORMAT(3x,13(i3,12x)) 
    204  9202 FORMAT(i3,41i3) 
    205  
    206  
    207       ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator) 
    208       ! ================= 
    209       ! define ahm1 and ahm2 at the right grid point position 
    210       ! (USER: modify ahm1 and ahm2 following your desiderata) 
    211        
    212        
    213       ! Decrease ahm to zahmeq m2/s in the tropics 
    214       ! (from 90 to 20 degre: ahm = constant 
    215       ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2 
    216       ! from  2.5 to  0 degre: ahm = constant 
    217       ! symmetric in the south hemisphere) 
    218  
    219       zahmeq = aht0 
    220        
    221       DO jj = 1, jpj 
    222          DO ji = 1, jpi 
    223             IF( ABS( gphif(ji,jj) ) >= 20. ) THEN 
    224                ahm2(ji,jj) =  ahm0 
    225             ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN 
    226                ahm2(ji,jj) =  zahmeq 
    227             ELSE 
    228                ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   & 
    229                   * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
    230             ENDIF 
    231             IF( ABS( gphit(ji,jj) ) >= 20. ) THEN 
    232                ahm1(ji,jj) =  ahm0 
    233             ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN 
    234                ahm1(ji,jj) =  zahmeq 
    235             ELSE 
    236                ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   & 
    237                   * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
    238             ENDIF 
    239          END DO 
    240       END DO 
    241  
    242       ! increase along western boundaries of equatorial strip 
    243       ! t-point 
    244       DO jj = 1, jpjm1 
    245          DO ji = 1, jpim1 
    246             zcoft = FLOAT( icof(ji,jj) ) / 100. 
    247             ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)  
    248          END DO 
    249       END DO 
    250       ! f-point 
    251       icof(:,:) = icof(:,:) * tmask(:,:,1) 
    252       DO jj = 1, jpjm1 
    253          DO ji = 1, jpim1   ! NO vector opt. 
    254             zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1) 
    255             IF( zmsk == 0. ) THEN 
    256                zcoff = 1. 
    257             ELSE 
    258                zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   & 
     238         ifreq = 40 
     239         il1 = 1 
     240         DO jn = 1, jpidta/ifreq+1 
     241            READ(inum,'(/)') 
     242            il2 = MIN( jpidta, il1+ifreq-1 ) 
     243            READ(inum,9201) ( ii, ji = il1, il2, 5 ) 
     244            READ(inum,'(/)') 
     245            DO jj = jpjdta, 1, -1 
     246               READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 ) 
     247            END DO 
     248            il1 = il1 + ifreq 
     249         END DO 
     250 
     251         DO jj = 1, nlcj 
     252            DO ji = 1, nlci 
     253               icof(ji,jj) = idata( mig(ji), mjg(jj) ) 
     254            END DO 
     255         END DO 
     256         DO jj = nlcj+1, jpj 
     257            DO ji = 1, nlci 
     258               icof(ji,jj) = icof(ji,nlcj) 
     259            END DO 
     260         END DO 
     261         DO jj = 1, jpj 
     262            DO ji = nlci+1, jpi 
     263               icof(ji,jj) = icof(nlci,jj) 
     264            END DO 
     265         END DO 
     266 
     2679101     FORMAT(1x,a15,2i8) 
     2689201     FORMAT(3x,13(i3,12x)) 
     2699202     FORMAT(i3,41i3) 
     270 
     271 
     272         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator) 
     273         ! ================= 
     274         ! define ahm1 and ahm2 at the right grid point position 
     275         ! (USER: modify ahm1 and ahm2 following your desiderata) 
     276 
     277 
     278         ! Decrease ahm to zahmeq m2/s in the tropics 
     279         ! (from 90 to 20 degre: ahm = constant 
     280         ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2 
     281         ! from  2.5 to  0 degre: ahm = constant 
     282         ! symmetric in the south hemisphere) 
     283 
     284         zahmeq = aht0 
     285 
     286         DO jj = 1, jpj 
     287            DO ji = 1, jpi 
     288               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN 
     289                  ahm2(ji,jj) =  ahm0 
     290               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN 
     291                  ahm2(ji,jj) =  zahmeq 
     292               ELSE 
     293                  ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   & 
     294                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
     295               ENDIF 
     296               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN 
     297                  ahm1(ji,jj) =  ahm0 
     298               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN 
     299                  ahm1(ji,jj) =  zahmeq 
     300               ELSE 
     301                  ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   & 
     302                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
     303               ENDIF 
     304            END DO 
     305         END DO 
     306 
     307         ! increase along western boundaries of equatorial strip 
     308         ! t-point 
     309         DO jj = 1, jpjm1 
     310            DO ji = 1, jpim1 
     311               zcoft = FLOAT( icof(ji,jj) ) / 100. 
     312               ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)  
     313            END DO 
     314         END DO 
     315         ! f-point 
     316         icof(:,:) = icof(:,:) * tmask(:,:,1) 
     317         DO jj = 1, jpjm1 
     318            DO ji = 1, jpim1   ! NO vector opt. 
     319               zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1) 
     320               IF( zmsk == 0. ) THEN 
     321                  zcoff = 1. 
     322               ELSE 
     323                  zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   & 
    259324                     / (zmsk * 100.) 
    260             ENDIF 
    261             ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj) 
    262          END DO 
    263       END DO 
    264 #endif 
     325               ENDIF 
     326               ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj) 
     327            END DO 
     328         END DO 
     329      ENDIF 
    265330       
    266331      ! Lateral boundary conditions on ( ahm1, ahm2 ) 
     
    323388      IF(lwp) WRITE(numout,*) '        orca_r1 configuration' 
    324389 
    325 #if defined key_antarctic 
    326 #     include "ldfdyn_antarctic.h90" 
    327 #elif defined key_arctic 
    328 #     include "ldfdyn_arctic.h90" 
    329 #else 
    330       ! Read 2d integer array to specify western boundary increase in the 
    331       ! ===================== equatorial strip (20N-20S) defined at t-points 
    332  
    333       CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   & 
    334          &           1, numout, lwp ) 
    335       REWIND inum 
    336       READ(inum,9101) clexp, iim, ijm 
    337       READ(inum,'(/)') 
    338       ifreq = 40 
    339       il1 = 1 
    340       DO jn = 1, jpidta/ifreq+1 
     390      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN 
     391! 
     392! 1.2 Modify ahm 
     393! -------------- 
     394         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean' 
     395         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm' 
     396         IF(lwp)WRITE(numout,*) '         north boundary increase' 
     397 
     398         ahm1(:,:) = ahm0 
     399         ahm2(:,:) = ahm0 
     400 
     401         ijpt0=max(1,min(49 -njmpp+1,jpj)) 
     402         ijpt1=max(0,min(49-njmpp+1,jpj-1)) 
     403         DO jj=ijpt0,ijpt1 
     404            ahm2(:,jj)=ahm0*2. 
     405            ahm1(:,jj)=ahm0*2. 
     406         END DO 
     407         ijpt0=max(1,min(48 -njmpp+1,jpj)) 
     408         ijpt1=max(0,min(48-njmpp+1,jpj-1)) 
     409         DO jj=ijpt0,ijpt1 
     410            ahm2(:,jj)=ahm0*1.9 
     411            ahm1(:,jj)=ahm0*1.75 
     412         END DO 
     413         ijpt0=max(1,min(47 -njmpp+1,jpj)) 
     414         ijpt1=max(0,min(47-njmpp+1,jpj-1)) 
     415         DO jj=ijpt0,ijpt1 
     416            ahm2(:,jj)=ahm0*1.5 
     417            ahm1(:,jj)=ahm0*1.25 
     418         END DO 
     419         ijpt0=max(1,min(46 -njmpp+1,jpj)) 
     420         ijpt1=max(0,min(46-njmpp+1,jpj-1)) 
     421         DO jj=ijpt0,ijpt1 
     422            ahm2(:,jj)=ahm0*1.1 
     423         END DO 
     424 
     425      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN 
     426! 1.2 Modify ahm  
     427! -------------- 
     428         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean' 
     429         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm' 
     430         IF(lwp)WRITE(numout,*) '         south and west boundary increase' 
     431 
     432 
     433         ahm1(:,:) = ahm0 
     434         ahm2(:,:) = ahm0 
     435 
     436         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj)) 
     437         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1)) 
     438         DO jj=ijpt0,ijpt1 
     439            ahm2(:,jj)=ahm0*2. 
     440            ahm1(:,jj)=ahm0*2. 
     441         END DO 
     442         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj)) 
     443         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1)) 
     444         DO jj=ijpt0,ijpt1 
     445            ahm2(:,jj)=ahm0*1.9 
     446            ahm1(:,jj)=ahm0*1.75 
     447         END DO 
     448         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj)) 
     449         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1)) 
     450         DO jj=ijpt0,ijpt1 
     451            ahm2(:,jj)=ahm0*1.5 
     452            ahm1(:,jj)=ahm0*1.25 
     453         END DO 
     454         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj)) 
     455         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1)) 
     456         DO jj=ijpt0,ijpt1 
     457            ahm2(:,jj)=ahm0*1.1 
     458         END DO 
     459      ELSE 
     460          
     461         ! Read 2d integer array to specify western boundary increase in the 
     462         ! ===================== equatorial strip (20N-20S) defined at t-points 
     463          
     464         CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   & 
     465            &           1, numout, lwp ) 
     466         REWIND inum 
     467         READ(inum,9101) clexp, iim, ijm 
    341468         READ(inum,'(/)') 
    342          il2 = MIN( jpidta, il1+ifreq-1 ) 
    343          READ(inum,9201) ( ii, ji = il1, il2, 5 ) 
    344          READ(inum,'(/)') 
    345          DO jj = jpjdta, 1, -1 
    346             READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 ) 
    347          END DO 
    348          il1 = il1 + ifreq 
    349       END DO 
    350        
    351       DO jj = 1, nlcj 
    352          DO ji = 1, nlci 
    353             icof(ji,jj) = idata( mig(ji), mjg(jj) ) 
    354          END DO 
    355       END DO 
    356       DO jj = nlcj+1, jpj 
    357          DO ji = 1, nlci 
    358             icof(ji,jj) = icof(ji,nlcj) 
    359          END DO 
    360       END DO 
    361       DO jj = 1, jpj 
    362          DO ji = nlci+1, jpi 
    363             icof(ji,jj) = icof(nlci,jj) 
    364          END DO 
    365       END DO 
    366  
    367  9101 FORMAT(1x,a15,2i8) 
    368  9201 FORMAT(3x,13(i3,12x)) 
    369  9202 FORMAT(i3,41i3) 
    370  
    371  
    372       ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator) 
    373       ! ================= 
    374       ! define ahm1 and ahm2 at the right grid point position 
    375       ! (USER: modify ahm1 and ahm2 following your desiderata) 
    376        
    377        
    378       ! Decrease ahm to zahmeq m2/s in the tropics 
    379       ! (from 90   to 20   degrees: ahm = scaled by local metrics 
    380       !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2 
    381       !  from  2.5 to  0   degrees: ahm = constant 
    382       ! symmetric in the south hemisphere) 
    383  
    384       zahmeq = aht0 
    385       zam20s = ahm0*COS( rad * 20. ) 
    386        
    387       DO jj = 1, jpj 
    388          DO ji = 1, jpi 
    389             IF( ABS( gphif(ji,jj) ) >= 20. ) THEN 
    390 !              leave as set in ldf_dyn_c2d 
    391             ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN 
    392                ahm2(ji,jj) =  zahmeq 
    393             ELSE 
    394                ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   & 
    395                   * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
    396             ENDIF 
    397             IF( ABS( gphit(ji,jj) ) >= 20. ) THEN 
    398 !             leave as set in ldf_dyn_c2d 
    399             ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN 
    400                ahm1(ji,jj) =  zahmeq 
    401             ELSE 
    402                ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   & 
    403                   * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
    404             ENDIF 
    405          END DO 
    406       END DO 
    407  
    408       ! increase along western boundaries of equatorial strip 
    409       ! t-point 
    410       DO jj = 1, jpjm1 
    411          DO ji = 1, jpim1 
    412           IF( ABS( gphit(ji,jj) ) < 20. ) THEN 
    413             zcoft = FLOAT( icof(ji,jj) ) / 100. 
    414             ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)  
    415           ENDIF 
    416          END DO 
    417       END DO 
    418       ! f-point 
    419       icof(:,:) = icof(:,:) * tmask(:,:,1) 
    420       DO jj = 1, jpjm1 
    421          DO ji = 1, jpim1 
    422           IF( ABS( gphif(ji,jj) ) < 20. ) THEN 
    423             zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1) 
    424             IF( zmsk == 0. ) THEN 
    425                zcoff = 1. 
    426             ELSE 
    427                zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   & 
    428                      / (zmsk * 100.) 
    429             ENDIF 
    430             ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj) 
    431           ENDIF 
    432          END DO 
    433       END DO 
    434 #endif 
     469         ifreq = 40 
     470         il1 = 1 
     471         DO jn = 1, jpidta/ifreq+1 
     472            READ(inum,'(/)') 
     473            il2 = MIN( jpidta, il1+ifreq-1 ) 
     474            READ(inum,9201) ( ii, ji = il1, il2, 5 ) 
     475            READ(inum,'(/)') 
     476            DO jj = jpjdta, 1, -1 
     477               READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 ) 
     478            END DO 
     479            il1 = il1 + ifreq 
     480         END DO 
     481 
     482         DO jj = 1, nlcj 
     483            DO ji = 1, nlci 
     484               icof(ji,jj) = idata( mig(ji), mjg(jj) ) 
     485            END DO 
     486         END DO 
     487         DO jj = nlcj+1, jpj 
     488            DO ji = 1, nlci 
     489               icof(ji,jj) = icof(ji,nlcj) 
     490            END DO 
     491         END DO 
     492         DO jj = 1, jpj 
     493            DO ji = nlci+1, jpi 
     494               icof(ji,jj) = icof(nlci,jj) 
     495            END DO 
     496         END DO 
     497 
     4989101     FORMAT(1x,a15,2i8) 
     4999201     FORMAT(3x,13(i3,12x)) 
     5009202     FORMAT(i3,41i3) 
     501 
     502 
     503         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator) 
     504         ! ================= 
     505         ! define ahm1 and ahm2 at the right grid point position 
     506         ! (USER: modify ahm1 and ahm2 following your desiderata) 
     507 
     508 
     509         ! Decrease ahm to zahmeq m2/s in the tropics 
     510         ! (from 90   to 20   degrees: ahm = scaled by local metrics 
     511         !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2 
     512         !  from  2.5 to  0   degrees: ahm = constant 
     513         ! symmetric in the south hemisphere) 
     514 
     515         zahmeq = aht0 
     516         zam20s = ahm0*COS( rad * 20. ) 
     517 
     518         DO jj = 1, jpj 
     519            DO ji = 1, jpi 
     520               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN 
     521                  !              leave as set in ldf_dyn_c2d 
     522               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN 
     523                  ahm2(ji,jj) =  zahmeq 
     524               ELSE 
     525                  ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   & 
     526                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
     527               ENDIF 
     528               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN 
     529                  !             leave as set in ldf_dyn_c2d 
     530               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN 
     531                  ahm1(ji,jj) =  zahmeq 
     532               ELSE 
     533                  ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   & 
     534                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
     535               ENDIF 
     536            END DO 
     537         END DO 
     538 
     539         ! increase along western boundaries of equatorial strip 
     540         ! t-point 
     541         DO jj = 1, jpjm1 
     542            DO ji = 1, jpim1 
     543               IF( ABS( gphit(ji,jj) ) < 20. ) THEN 
     544                  zcoft = FLOAT( icof(ji,jj) ) / 100. 
     545                  ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)  
     546               ENDIF 
     547            END DO 
     548         END DO 
     549         ! f-point 
     550         icof(:,:) = icof(:,:) * tmask(:,:,1) 
     551         DO jj = 1, jpjm1 
     552            DO ji = 1, jpim1 
     553               IF( ABS( gphif(ji,jj) ) < 20. ) THEN 
     554                  zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1) 
     555                  IF( zmsk == 0. ) THEN 
     556                     zcoff = 1. 
     557                  ELSE 
     558                     zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   & 
     559                        / (zmsk * 100.) 
     560                  ENDIF 
     561                  ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj) 
     562               ENDIF 
     563            END DO 
     564         END DO 
     565      ENDIF 
    435566       
    436567      ! Lateral boundary conditions on ( ahm1, ahm2 ) 
Note: See TracChangeset for help on using the changeset viewer.