Ignore:
Timestamp:
Nov 21, 2019, 4:43:45 PM (5 years ago)
Author:
lguez
Message:

Merge revisions 3427:3600 of trunk into branch Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Ocean_skin

  • LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90

    r3102 r3605  
    9797  !$OMP THREADPRIVATE(cpl_atm_co22D)
    9898
     99!!!!!!!!!! variable for calving
     100  INTEGER, PARAMETER :: nb_zone_calving = 3
     101  REAL,ALLOCATABLE, DIMENSION(:,:,:),SAVE :: area_calving
     102  !$OMP THREADPRIVATE(area_calving)
     103  REAL,ALLOCATABLE, DIMENSION(:,:),SAVE :: cell_area2D
     104  !$OMP THREADPRIVATE(cell_area2D)
     105  INTEGER, SAVE :: ind_calving(nb_zone_calving)
     106  !$OMP THREADPRIVATE(ind_calving)
     107
     108  LOGICAL,SAVE :: cpl_old_calving
     109  !$OMP THREADPRIVATE(cpl_old_calving)
     110 
    99111CONTAINS
    100112!
    … …  
    105117    USE surface_data
    106118    USE indice_sol_mod
    107     USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
     119    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat, grid1dTo2d_glo, klon_glo, grid_type, unstructured, regular_lonlat
    108120    USE time_phylmdz_mod, ONLY: annee_ref, day_ini, itau_phy, itaufin_phy
    109121    USE print_control_mod, ONLY: lunout
     122    USE geometry_mod, ONLY : longitude_deg, latitude_deg, ind_cell_glo, cell_area
     123    USE ioipsl_getin_p_mod, ONLY: getin_p
    110124
    111125! Input arguments
    … …  
    127141    CHARACTER(len = 80)               :: abort_message
    128142    CHARACTER(len=80)                 :: clintocplnam, clfromcplnam
     143    REAL, DIMENSION(klon_mpi)         :: rlon_mpi, rlat_mpi, cell_area_mpi
     144    INTEGER, DIMENSION(klon_mpi)           :: ind_cell_glo_mpi
     145    REAL, DIMENSION(nbp_lon,jj_nb)         :: lon2D, lat2D
     146    INTEGER :: mask_calving(nbp_lon,jj_nb,nb_zone_calving)
     147    REAL :: pos
     148
     149!***************************************
     150! Use old calving or not (default new calving method)
     151! New calving method should be used with DYNAMICO and when using new coupling
     152! weights.
     153    cpl_old_calving=.FALSE.
     154    CALL getin_p("cpl_old_calving",cpl_old_calving)
     155
    129156
    130157!*************************************************************************************
    … …  
    204231
    205232! Allocate variable in carbon_cycle_mod
    206        ALLOCATE(fco2_ocn_day(klon), stat = error)
     233       IF (.NOT.ALLOCATED(fco2_ocn_day)) ALLOCATE(fco2_ocn_day(klon), stat = error)
    207234       sum_error = sum_error + error
    208     END IF
    209 
     235    ENDIF
     236
     237! calving initialization
     238    ALLOCATE(area_calving(nbp_lon, jj_nb, nb_zone_calving), stat = error)
     239    sum_error = sum_error + error
     240    ALLOCATE(cell_area2D(nbp_lon, jj_nb), stat = error)   
     241    sum_error = sum_error + error
     242
     243
     244    CALL gather_omp(longitude_deg,rlon_mpi)
     245    CALL gather_omp(latitude_deg,rlat_mpi)
     246    CALL gather_omp(ind_cell_glo,ind_cell_glo_mpi)
     247    CALL gather_omp(cell_area,cell_area_mpi)
     248     
     249    IF (is_omp_master) THEN
     250      CALL Grid1DTo2D_mpi(rlon_mpi,lon2D)
     251      CALL Grid1DTo2D_mpi(rlat_mpi,lat2D)
     252      CALL Grid1DTo2D_mpi(cell_area_mpi,cell_area2D)
     253      mask_calving(:,:,:) = 0
     254      WHERE ( lat2D >= 40) mask_calving(:,:,1) = 1
     255      WHERE ( lat2D < 40 .AND. lat2D > -50) mask_calving(:,:,2) = 1
     256      WHERE ( lat2D <= -50) mask_calving(:,:,3) = 1
     257   
     258   
     259      DO i=1,nb_zone_calving
     260        area_calving(:,:,i)=mask_calving(:,:,i)*cell_area2D(:,:)
     261        pos=1
     262        IF (i>1) pos = 1 + ((nbp_lon*nbp_lat-1)*(i-1))/(nb_zone_calving-1)
     263     
     264        ind_calving(i)=0
     265        IF (grid_type==unstructured) THEN
     266
     267          DO ig=1,klon_mpi
     268            IF (ind_cell_glo_mpi(ig)==pos) ind_calving(i)=ig
     269          ENDDO
     270
     271        ELSE IF (grid_type==regular_lonlat) THEN
     272          IF ((ij_begin<=pos .AND. ij_end>=pos) .OR. (ij_begin<=pos .AND. is_south_pole_dyn )) THEN
     273            ind_calving(i)=pos-(jj_begin-1)*nbp_lon
     274          ENDIF
     275        ENDIF
     276     
     277      ENDDO
     278    ENDIF
     279   
     280           
    210281    IF (sum_error /= 0) THEN
    211282       abort_message='Pb allocation variables couplees'
    … …  
    236307       idayref = day_ini
    237308       CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
    238        CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlon,zx_lon)
     309       CALL grid1dTo2d_glo(rlon,zx_lon)
    239310       DO i = 1, nbp_lon
    240311          zx_lon(i,1) = rlon(i+1)
    241312          zx_lon(i,nbp_lat) = rlon(i+1)
    242313       ENDDO
    243        CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlat,zx_lat)
     314       CALL grid1dTo2d_glo(rlat,zx_lat)
    244315       clintocplnam="cpl_atm_tauflx"
    245316       CALL histbeg(clintocplnam,nbp_lon,zx_lon(:,1),nbp_lat,zx_lat(1,:),&
    … …  
    259330                "-",nbp_lon,nbp_lat,nhoridct,1,1,1,-99,32,"inst",dtime,dtime)
    260331         ENDIF
    261        END DO
     332       ENDDO
    262333       CALL histend(nidct)
    263334       CALL histsync(nidct)
    … …  
    272343                "-",nbp_lon,nbp_lat,nhoridcs,1,1,1,-99,32,"inst",dtime,dtime)
    273344         ENDIF
    274        END DO
     345       ENDDO
    275346       CALL histend(nidcs)
    276347       CALL histsync(nidcs)
    … …  
    286357       abort_message='carbon_cycle_cpl does not work with opa8'
    287358       CALL abort_physic(modname,abort_message,1)
    288     END IF
     359    ENDIF
    289360
    290361  END SUBROUTINE cpl_init
    … …  
    356427                CALL histwrite(nidcs,inforecv(i)%name,itau_w,tab_read_flds(:,:,i),nbp_lon*(nbp_lat),ndexcs)
    357428            ENDIF
    358           END DO
     429          ENDDO
    359430       ENDIF
    360431
    … …  
    415486       ENDDO
    416487
    417     END IF ! if time to receive
     488    ENDIF ! if time to receive
    418489
    419490  END SUBROUTINE cpl_receive_frac
    … …  
    466537       DO i=1,klon
    467538          index(i)=i
    468        END DO
     539       ENDDO
    469540       CALL cpl2gath(read_co2, fco2_ocn_day, klon, index)
    470     END IF
     541    ENDIF
    471542
    472543!*************************************************************************************
    … …  
    477548    DO i=1, knon
    478549       tsurf_new(i) = tsurf_new(i)/(1. - sic_new(i))
    479     END DO
     550    ENDDO
    480551
    481552  END SUBROUTINE cpl_receive_ocean_fields
    … …  
    529600       tsurf_new(i) = tsurf_new(i) / sic_new(i)
    530601       alb_new(i)   = alb_new(i)   / sic_new(i)
    531     END DO
     602    ENDDO
    532603
    533604  END SUBROUTINE cpl_receive_seaice_fields
    … …  
    637708          cpl_atm_co2(ig,cpl_index) = cpl_atm_co2(ig,cpl_index) + &
    638709               co2_send(knindex(ig))/ REAL(nexca)
    639        END IF
     710!!---OB: this is correct but why knindex ??
     711       ENDIF
    640712     ENDDO
    641713
    … …  
    682754             ALLOCATE(cpl_atm_co22D(nbp_lon,jj_nb), stat=error)
    683755             sum_error = sum_error + error
    684           END IF
     756          ENDIF
    685757
    686758          IF (sum_error /= 0) THEN
    … …  
    886958             ALLOCATE(cpl_atm_co22D(nbp_lon,jj_nb), stat=error)
    887959             sum_error = sum_error + error
    888           END IF
     960          ENDIF
    889961
    890962          IF (sum_error /= 0) THEN
    … …  
    917989       DO ig = 1, knon
    918990          cpl_fder_tmp(knindex(ig))=cpl_fder(ig,cpl_index)
    919        END DO
     991       ENDDO
    920992       CALL gath2cpl(cpl_fder_tmp(:), cpl_fder2D(:,:,cpl_index), &
    921993            klon, unity)
    … …  
    10851157! Local variables
    10861158!*************************************************************************************
    1087     INTEGER                                              :: error, sum_error, j
     1159    INTEGER                                              :: error, sum_error, i,j,k
    10881160    INTEGER                                              :: itau_w
    10891161    INTEGER                                              :: time_sec
    … …  
    11021174! Table with all fields to send to coupler
    11031175    REAL, DIMENSION(nbp_lon, jj_nb, maxsend)                 :: tab_flds
    1104     REAL, DIMENSION(klon_mpi)                            :: rlon_mpi, rlat_mpi
    1105 
     1176    REAL, DIMENSION(klon_mpi)                                :: rlon_mpi, rlat_mpi
     1177    REAL  :: calving(nb_zone_calving)
     1178    REAL  :: calving_glo(nb_zone_calving)
     1179   
    11061180#ifdef CPP_MPI
    11071181    INCLUDE 'mpif.h'
    … …  
    11301204   
    11311205    IF (version_ocean=='nemo') THEN
    1132        tab_flds(:,:,ids_liqrun) = cpl_rriv2D(:,:) + cpl_rcoa2D(:,:)
     1206       tab_flds(:,:,ids_liqrun) = (cpl_rriv2D(:,:) + cpl_rcoa2D(:,:))
    11331207       IF (carbon_cycle_cpl) tab_flds(:,:,ids_atmco2)=cpl_atm_co22D(:,:)
    11341208    ELSE IF (version_ocean=='opa8') THEN
    … …  
    11391213       tab_flds(:,:,ids_runcoa) = cpl_rcoa2D(:,:)
    11401214       tab_flds(:,:,ids_rivflu) = cpl_rriv2D(:,:)
    1141     END IF
     1215    ENDIF
    11421216
    11431217!*************************************************************************************
    … …  
    11581232    IF (is_omp_root) THEN
    11591233
    1160       DO j = 1, jj_nb
    1161          tmp_calv(:,j) = DOT_PRODUCT (cpl_rlic2D(1:nbp_lon,j), &
    1162               pctsrf2D(1:nbp_lon,j,is_lic)) / REAL(nbp_lon)
    1163       ENDDO
    1164    
    1165    
    1166       IF (is_parallel) THEN
    1167          IF (.NOT. is_north_pole_dyn) THEN
     1234      IF (cpl_old_calving) THEN   ! use old calving
     1235
     1236        DO j = 1, jj_nb
     1237           tmp_calv(:,j) = DOT_PRODUCT (cpl_rlic2D(1:nbp_lon,j), &
     1238                pctsrf2D(1:nbp_lon,j,is_lic)) / REAL(nbp_lon)
     1239        ENDDO
     1240   
     1241   
     1242        IF (is_parallel) THEN
     1243           IF (.NOT. is_north_pole_dyn) THEN
    11681244#ifdef CPP_MPI
    1169             CALL MPI_RECV(Up,1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,status,error)
    1170             CALL MPI_SEND(tmp_calv(1,1),1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,error)
     1245              CALL MPI_RECV(Up,1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,status,error)
     1246              CALL MPI_SEND(tmp_calv(1,1),1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,error)
    11711247#endif
    1172          ENDIF
     1248           ENDIF
    11731249       
    1174          IF (.NOT. is_south_pole_dyn) THEN
     1250           IF (.NOT. is_south_pole_dyn) THEN
    11751251#ifdef CPP_MPI
    1176             CALL MPI_SEND(tmp_calv(1,jj_nb),1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,error)
    1177             CALL MPI_RECV(down,1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,status,error)
     1252              CALL MPI_SEND(tmp_calv(1,jj_nb),1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,error)
     1253              CALL MPI_RECV(down,1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,status,error)
    11781254#endif
    1179          ENDIF
     1255           ENDIF
    11801256         
    1181          IF (.NOT. is_north_pole_dyn .AND. ii_begin /=1) THEN
    1182             Up=Up+tmp_calv(nbp_lon,1)
    1183             tmp_calv(:,1)=Up
    1184          ENDIF
     1257           IF (.NOT. is_north_pole_dyn .AND. ii_begin /=1) THEN
     1258              Up=Up+tmp_calv(nbp_lon,1)
     1259              tmp_calv(:,1)=Up
     1260           ENDIF
     1261           
     1262           IF (.NOT. is_south_pole_dyn .AND. ii_end /= nbp_lon) THEN
     1263              Down=Down+tmp_calv(1,jj_nb)
     1264              tmp_calv(:,jj_nb)=Down
     1265           ENDIF
     1266        ENDIF
     1267        tab_flds(:,:,ids_calvin) = tmp_calv(:,:)
     1268
     1269      ELSE
     1270         ! cpl_old_calving=FALSE
     1271         ! To be used with new method for calculation of coupling weights
     1272         DO k=1,nb_zone_calving
     1273            calving(k)=0
     1274            DO j = 1, jj_nb
     1275               calving(k)= calving(k)+DOT_PRODUCT(cpl_rlic2D(:,j)*area_calving(:,j,k),pctsrf2D(:,j,is_lic))
     1276            ENDDO
     1277         ENDDO
    11851278         
    1186          IF (.NOT. is_south_pole_dyn .AND. ii_end /= nbp_lon) THEN
    1187             Down=Down+tmp_calv(1,jj_nb)
    1188             tmp_calv(:,jj_nb)=Down       
    1189          ENDIF
     1279#ifdef CPP_MPI
     1280         CALL MPI_ALLREDUCE(calving, calving_glo, nb_zone_calving, MPI_REAL_LMDZ, MPI_SUM, COMM_LMDZ_PHY, error)
     1281#endif
     1282         
     1283         tab_flds(:,:,ids_calvin) = 0
     1284         DO k=1,nb_zone_calving
     1285            IF (ind_calving(k)>0 ) THEN
     1286               j=(ind_calving(k)-1)/nbp_lon + 1
     1287               i=MOD(ind_calving(k)-1,nbp_lon)+1
     1288               tab_flds(i,j,ids_calvin) = calving_glo(k)
     1289            ENDIF
     1290         ENDDO
     1291         
    11901292      ENDIF
    11911293     
    1192       tab_flds(:,:,ids_calvin) = tmp_calv(:,:)
    1193 
    11941294!*************************************************************************************
    11951295! Calculate total flux for snow, rain and wind with weighted addition using the
    … …  
    12521352                  cpl_tauy2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
    12531353          ENDWHERE
    1254        END IF
     1354       ENDIF
    12551355
    12561356    ENDIF ! is_omp_root
    … …  
    13361436       DEALLOCATE(cpl_atm_co22D, stat=error )
    13371437       sum_error = sum_error + error
    1338     END IF
     1438    ENDIF
    13391439
    13401440    IF (sum_error /= 0) THEN
Note: See TracChangeset for help on using the changeset viewer.