source: LMDZ6/trunk/libf/phylmd/phyaqua_mod.F90 @ 3531

Last change on this file since 3531 was 3531, checked in by Laurent Fairhead, 5 years ago

Replaced STOP statements by a call to abort_physic in phylmd as per ticket #86
Still some work to be done in phylmd subdirectories
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 30.9 KB
Line 
1MODULE phyaqua_mod
2  ! Routines complementaires pour la physique planetaire.
3  IMPLICIT NONE
4
5CONTAINS
6
7  SUBROUTINE iniaqua(nlon, iflag_phys)
8
9    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10    ! Creation d'un etat initial et de conditions aux limites
11    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
12    ! du modele LMDZ dans sa version terrestre.
13    ! iflag_phys est un parametre qui controle
14    ! iflag_phys = N
15    ! de 100 a 199 : aqua planetes avec SST forcees
16    ! N-100 determine le type de SSTs
17    ! de 200 a 299 : terra planetes avec Ts calcule
18
19    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20
21    USE dimphy, ONLY: klon
22    USE geometry_mod, ONLY : latitude
23    USE surface_data, ONLY: type_ocean, ok_veget
24    USE pbl_surface_mod, ONLY: pbl_surface_init
25    USE fonte_neige_mod, ONLY: fonte_neige_init
26    USE phys_state_var_mod
27    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
28                                day_ini,day_end
29    USE indice_sol_mod
30    USE nrtype, ONLY: pi
31!    USE ioipsl
32    USE mod_phys_lmdz_para, ONLY: is_master
33    USE mod_phys_lmdz_transfert_para, ONLY: bcast
34    USE mod_grid_phy_lmdz
35    USE ioipsl_getin_p_mod, ONLY : getin_p
36    IMPLICIT NONE
37
38    include "YOMCST.h"
39    include "clesphys.h"
40    include "dimsoil.h"
41
42    INTEGER, INTENT (IN) :: nlon, iflag_phys
43    ! IM ajout latfi, lonfi
44!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
45
46    INTEGER type_profil, type_aqua
47
48    ! Ajouts initialisation des surfaces
49    REAL :: run_off_lic_0(nlon)
50    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
51    REAL :: tsoil(nlon, nsoilmx, nbsrf)
52    REAL :: tslab(nlon), seaice(nlon)
53    REAL fder(nlon)
54
55
56
57    ! Arguments :
58    ! -----------
59
60    ! integer radpas
61    INTEGER it, unit, i, k, itap
62
63    REAL rugos, albedo
64    REAL tsurf
65    REAL time, timestep, day, day0
66    REAL qsol_f
67    REAL rugsrel(nlon)
68    LOGICAL alb_ocean
69
70    CHARACTER *80 ans, file_forctl, file_fordat, file_start
71    CHARACTER *100 file, var
72    CHARACTER *2 cnbl
73
74    REAL phy_nat(nlon, 360)
75    REAL phy_alb(nlon, 360)
76    REAL phy_sst(nlon, 360)
77    REAL phy_bil(nlon, 360)
78    REAL phy_rug(nlon, 360)
79    REAL phy_ice(nlon, 360)
80    REAL phy_fter(nlon, 360)
81    REAL phy_foce(nlon, 360)
82    REAL phy_fsic(nlon, 360)
83    REAL phy_flic(nlon, 360)
84
85    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
86!$OMP THREADPRIVATE(read_climoz)
87
88    ! -------------------------------------------------------------------------
89    ! declaration pour l'appel a phyredem
90    ! -------------------------------------------------------------------------
91
92    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
93    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
94    ! real pbl_tke(nlon,llm,nbsrf)
95    ! real rain_fall(nlon),snow_fall(nlon)
96    ! real solsw(nlon), sollw(nlon),radsol(nlon)
97    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
98    ! real ratqs(nlon,llm)
99    ! real clwcon(nlon,llm)
100
101    INTEGER longcles
102    PARAMETER (longcles=20)
103    REAL clesphy0(longcles)
104
105
106    ! -----------------------------------------------------------------------
107    ! dynamial tendencies :
108    ! ---------------------
109
110    INTEGER l, ierr, aslun
111
112    REAL paire
113
114    ! Local
115    CHARACTER (LEN=20) :: modname='phyaqua'
116    CHARACTER (LEN=80) :: abort_message
117
118
119    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
120    ! INITIALISATIONS
121    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
122
123    ! -----------------------------------------------------------------------
124    ! Initialisations  des constantes
125    ! -------------------------------
126
127
128    type_aqua = iflag_phys/100
129    type_profil = iflag_phys - type_aqua*100
130    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
131
132    IF (klon/=nlon) THEN
133      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
134      abort_message= 'probleme de dimensions dans iniaqua'
135      CALL abort_physic(modname,abort_message,1)
136    END IF
137    CALL phys_state_var_init(read_climoz)
138
139
140    read_climoz = 0
141    day0 = 217.
142    day = day0
143    it = 0
144    time = 0.
145
146    ! -----------------------------------------------------------------------
147    ! initialisations de la physique
148    ! -----------------------------------------------------------------------
149
150    day_ini = day_ref
151    day_end = day_ini + ndays
152
153    nbapp_rad = 24
154    CALL getin_p('nbapp_rad', nbapp_rad)
155
156    ! ---------------------------------------------------------------------
157    ! Creation des conditions aux limites:
158    ! ------------------------------------
159    ! Initialisations des constantes
160    ! Ajouter les manquants dans planete.def... (albedo etc)
161    co2_ppm = 348.
162    CALL getin_p('co2_ppm', co2_ppm)
163
164    solaire = 1365.
165    CALL getin_p('solaire', solaire)
166 
167    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
168    ! type_aqua
169    alb_ocean = .TRUE.
170    CALL getin_p('alb_ocean', alb_ocean)
171
172    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
173    WRITE (*, *) 'iniaqua: solaire=', solaire
174    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
175
176    radsol = 0.
177    qsol_f = 10.
178
179    ! Conditions aux limites:
180    ! -----------------------
181
182    qsol(:) = qsol_f
183    rugsrel = 0.0 ! (rugsrel = rugoro)
184    rugoro = 0.0
185    u_ancien = 0.0
186    v_ancien = 0.0
187    agesno = 50.0
188    ! Relief plat
189    zmea = 0.
190    zstd = 0.
191    zsig = 0.
192    zgam = 0.
193    zthe = 0.
194    zpic = 0.
195    zval = 0.
196
197    ! Une seule surface
198    pctsrf = 0.
199    IF (type_aqua==1) THEN
200      rugos = 1.E-4
201      albedo = 0.19
202      pctsrf(:, is_oce) = 1.
203    ELSE IF (type_aqua==2) THEN
204      rugos = 0.03
205      albedo = 0.1
206      pctsrf(:, is_ter) = 1.
207    END IF
208
209    CALL getin_p('rugos', rugos)
210
211    WRITE (*, *) 'iniaqua: rugos=', rugos
212    zmasq(:) = pctsrf(:, is_ter)
213
214    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
215    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
216
217    ! Si alb_ocean on calcule un albedo oceanique moyen
218    ! if (alb_ocean) then
219    ! Voir pourquoi on avait ca.
220    ! CALL ini_alb_oce(phy_alb)
221    ! else
222    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
223    ! endif !alb_ocean
224
225    DO i = 1, 360
226      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
227      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
228
229      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
230      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
231      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
232      phy_ice(:, i) = 0.0 ! fraction de glace (?)
233      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
234      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
235      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
236      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
237    END DO
238    ! IM calcul profil sst
239    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
240
241    IF (grid_type==unstructured) THEN
242      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
243                             phy_fter, phy_foce, phy_flic, phy_fsic)
244    ELSE
245     
246       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
247                     phy_fter, phy_foce, phy_flic, phy_fsic)
248    ENDIF
249
250    ! ---------------------------------------------------------------------
251    ! Ecriture de l'etat initial:
252    ! ---------------------------
253
254
255    ! Ecriture etat initial physique
256
257    timestep = pdtphys
258    radpas = nint(rday/timestep/float(nbapp_rad))
259
260    DO i = 1, longcles
261      clesphy0(i) = 0.
262    END DO
263    clesphy0(1) = float(iflag_con)
264    clesphy0(2) = float(nbapp_rad)
265    ! IF( cycle_diurne  ) clesphy0(3) =  1.
266    clesphy0(3) = 1. ! cycle_diurne
267    clesphy0(4) = 1. ! soil_model
268    clesphy0(5) = 1. ! new_oliq
269    clesphy0(6) = 0. ! ok_orodr
270    clesphy0(7) = 0. ! ok_orolf
271    clesphy0(8) = 0. ! ok_limitvrai
272
273
274    ! =======================================================================
275    ! Profils initiaux
276    ! =======================================================================
277
278    ! On initialise les temperatures de surfaces comme les sst
279    DO i = 1, nlon
280      ftsol(i, :) = phy_sst(i, 1)
281      tsoil(i, :, :) = phy_sst(i, 1)
282      tslab(i) = phy_sst(i, 1)
283    END DO
284
285    falbe(:, :) = albedo
286    falblw(:, :) = albedo
287    rain_fall(:) = 0.
288    snow_fall(:) = 0.
289    solsw(:) = 0.
290    sollw(:) = 0.
291    radsol(:) = 0.
292
293    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
294    ! intialisation bidon mais pas grave
295    t_ancien(:, :) = 0.
296    q_ancien(:, :) = 0.
297    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
298    rnebcon = 0.
299    ratqs = 0.
300    clwcon = 0.
301    pbl_tke = 1.E-8
302
303    ! variables supplementaires pour appel a plb_surface_init
304    fder(:) = 0.
305    seaice(:) = 0.
306    run_off_lic_0 = 0.
307    fevap = 0.
308
309
310    ! Initialisations necessaires avant phyredem
311    type_ocean = 'force'
312    CALL fonte_neige_init(run_off_lic_0)
313    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
314    snsrf(:, :) = 0. ! couverture de neige des sous surface
315    z0m(:, :) = rugos ! couverture de neige des sous surface
316    z0h=z0m
317
318
319    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
320
321    PRINT *, 'iniaqua: before phyredem'
322
323    pbl_tke(:,:,:) = 1.e-8
324    falb1 = albedo
325    falb2 = albedo
326    zmax0 = 0.
327    f0 = 0.
328    sig1 = 0.
329    w01 = 0.
330    wake_deltat = 0.
331    wake_deltaq = 0.
332    wake_s = 0.
333    wake_dens = 0.
334    wake_cstar = 0.
335    wake_pe = 0.
336    wake_fip = 0.
337    fm_therm = 0.
338    entr_therm = 0.
339    detr_therm = 0.
340    ale_bl = 0.
341    ale_bl_trig =0.
342    alp_bl =0.
343    treedrg(:,:,:)=0.
344
345!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
346!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
347!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
348    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
349    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
350
351!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
352!ym probably the uninitialized value was 0 for standard (regular grid) case
353    falb_dif(:,:,:)=0
354
355
356    CALL phyredem('startphy.nc')
357
358    PRINT *, 'iniaqua: after phyredem'
359    CALL phys_state_var_end
360
361    RETURN
362  END SUBROUTINE iniaqua
363
364
365  ! ====================================================================
366  ! ====================================================================
367  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
368    USE dimphy
369    IMPLICIT NONE
370    ! ====================================================================
371    ! =============================================================
372    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
373    ! Auteur : A. Campoy et F. Hourdin
374    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
375    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
376    ! connaissant la declinaison, la latitude et la longitude.
377
378    ! Dans cette version particuliere, on calcule le rayonnement
379    ! moyen sur l'année à chaque latitude.
380    ! angle zenithal calculé pour obtenir un
381    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
382    ! en moyenne annuelle.
383    ! Spécifique de la terre. Utilisé pour les aqua planetes.
384
385    ! Rque   : Different de la routine angle en ce sens que zenang
386    ! fournit des moyennes de pmu0 et non des valeurs
387    ! instantanees, du coup frac prend toutes les valeurs
388    ! entre 0 et 1.
389    ! Date   : premiere version le 13 decembre 1994
390    ! revu pour  GCM  le 30 septembre 1996
391    ! ===============================================================
392    ! longi----INPUT : la longitude vraie de la terre dans son plan
393    ! solaire a partir de l'equinoxe de printemps (degre)
394    ! gmtime---INPUT : temps universel en fraction de jour
395    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
396    ! lat------INPUT : latitude en degres
397    ! long-----INPUT : longitude en degres
398    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
399    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
400    ! ================================================================
401    include "YOMCST.h"
402    ! ================================================================
403    LOGICAL cycle_diurne
404    REAL gmtime
405    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
406    ! ================================================================
407    INTEGER i
408    REAL gmtime1, gmtime2
409    REAL pi_local
410
411
412    REAL rmu0m(klon), rmu0a(klon)
413
414
415    pi_local = 4.0*atan(1.0)
416
417    ! ================================================================
418    ! Calcul de l'angle zenithal moyen sur la journee
419    ! ================================================================
420
421    DO i = 1, klon
422      fract(i) = 1.
423      ! Calcule du flux moyen
424      IF (abs(rlat(i))<=28.75) THEN
425        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
426      ELSE IF (abs(rlat(i))<=43.75) THEN
427        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
428      ELSE IF (abs(rlat(i))<=71.25) THEN
429        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
430      ELSE
431        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
432      END IF
433    END DO
434
435    ! ================================================================
436    ! Avec ou sans cycle diurne
437    ! ================================================================
438
439    IF (cycle_diurne) THEN
440
441      ! On redecompose flux  au sommet suivant un cycle diurne idealise
442      ! identique a toutes les latitudes.
443
444      DO i = 1, klon
445        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
446        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
447          rmu0a(i)/sqrt(2.)
448      END DO
449
450      DO i = 1, klon
451        IF (rmu0(i)<=0.) THEN
452          rmu0(i) = 0.
453          fract(i) = 0.
454        ELSE
455          fract(i) = 1.
456        END IF
457      END DO
458
459      ! Affichage de l'angel zenitale
460      ! print*,'************************************'
461      ! print*,'************************************'
462      ! print*,'************************************'
463      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
464      ! print*,'rmu0m=',rmu0m(i)
465      ! print*,'rmu0a=',rmu0a(i)
466      ! print*,'rmu0=',rmu0(i)
467
468    ELSE
469
470      DO i = 1, klon
471        fract(i) = 0.5
472        rmu0(i) = rmu0m(i)*2.
473      END DO
474
475    END IF
476
477    RETURN
478  END SUBROUTINE zenang_an
479
480  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
481
482  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
483      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
484
485    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
486    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
487#ifdef CPP_XIOS
488    USE xios
489#endif
490    IMPLICIT NONE
491
492    include "netcdf.inc"
493
494    INTEGER, INTENT (IN) :: klon
495    REAL, INTENT (IN) :: phy_nat(klon, 360)
496    REAL, INTENT (IN) :: phy_alb(klon, 360)
497    REAL, INTENT (IN) :: phy_sst(klon, 360)
498    REAL, INTENT (IN) :: phy_bil(klon, 360)
499    REAL, INTENT (IN) :: phy_rug(klon, 360)
500    REAL, INTENT (IN) :: phy_ice(klon, 360)
501    REAL, INTENT (IN) :: phy_fter(klon, 360)
502    REAL, INTENT (IN) :: phy_foce(klon, 360)
503    REAL, INTENT (IN) :: phy_flic(klon, 360)
504    REAL, INTENT (IN) :: phy_fsic(klon, 360)
505
506    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
507      ! on the whole physics grid
508 
509#ifdef CPP_XIOS
510    PRINT *, 'writelim: Ecriture du fichier limit'
511
512    CALL gather_omp(phy_foce, phy_mpi)
513    IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
514
515    CALL gather_omp(phy_fsic, phy_mpi)
516    IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
517     
518    CALL gather_omp(phy_fter, phy_mpi)
519    IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
520     
521    CALL gather_omp(phy_flic, phy_mpi)
522    IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
523
524    CALL gather_omp(phy_sst, phy_mpi)
525    IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
526
527    CALL gather_omp(phy_bil, phy_mpi)
528    IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
529
530    CALL gather_omp(phy_alb, phy_mpi)
531    IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
532
533    CALL gather_omp(phy_rug, phy_mpi)
534    IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
535#endif
536  END SUBROUTINE writelim_unstruct
537
538
539
540  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
541      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
542
543    USE mod_phys_lmdz_para, ONLY: is_master
544    USE mod_grid_phy_lmdz, ONLY: klon_glo
545    USE mod_phys_lmdz_transfert_para, ONLY: gather
546    IMPLICIT NONE
547    include "netcdf.inc"
548
549    INTEGER, INTENT (IN) :: klon
550    REAL, INTENT (IN) :: phy_nat(klon, 360)
551    REAL, INTENT (IN) :: phy_alb(klon, 360)
552    REAL, INTENT (IN) :: phy_sst(klon, 360)
553    REAL, INTENT (IN) :: phy_bil(klon, 360)
554    REAL, INTENT (IN) :: phy_rug(klon, 360)
555    REAL, INTENT (IN) :: phy_ice(klon, 360)
556    REAL, INTENT (IN) :: phy_fter(klon, 360)
557    REAL, INTENT (IN) :: phy_foce(klon, 360)
558    REAL, INTENT (IN) :: phy_flic(klon, 360)
559    REAL, INTENT (IN) :: phy_fsic(klon, 360)
560
561    REAL :: phy_glo(klon_glo, 360) ! temporary variable, to store phy_***(:)
562      ! on the whole physics grid
563    INTEGER :: k
564    INTEGER ierr
565    INTEGER dimfirst(3)
566    INTEGER dimlast(3)
567
568    INTEGER nid, ndim, ntim
569    INTEGER dims(2), debut(2), epais(2)
570    INTEGER id_tim
571    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
572    INTEGER id_fter, id_foce, id_fsic, id_flic
573
574    IF (is_master) THEN
575
576      PRINT *, 'writelim: Ecriture du fichier limit'
577
578      ierr = nf_create('limit.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
579
580      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
581        'Fichier conditions aux limites')
582      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
583      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
584      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
585
586      dims(1) = ndim
587      dims(2) = ntim
588
589#ifdef NC_DOUBLE
590      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
591#else
592      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
593#endif
594      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
595
596#ifdef NC_DOUBLE
597      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
598#else
599      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
600#endif
601      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
602        'Nature du sol (0,1,2,3)')
603
604#ifdef NC_DOUBLE
605      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
606#else
607      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
608#endif
609      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
610        'Temperature superficielle de la mer')
611
612#ifdef NC_DOUBLE
613      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
614#else
615      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
616#endif
617      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
618        'Reference flux de chaleur au sol')
619
620#ifdef NC_DOUBLE
621      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
622#else
623      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
624#endif
625      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
626
627#ifdef NC_DOUBLE
628      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
629#else
630      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
631#endif
632      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
633
634#ifdef NC_DOUBLE
635      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
636#else
637      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
638#endif
639      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
640#ifdef NC_DOUBLE
641      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
642#else
643      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
644#endif
645      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
646#ifdef NC_DOUBLE
647      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
648#else
649      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
650#endif
651      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
652#ifdef NC_DOUBLE
653      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
654#else
655      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
656#endif
657      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
658
659      ierr = nf_enddef(nid)
660      IF (ierr/=nf_noerr) THEN
661        WRITE (*, *) 'writelim error: failed to end define mode'
662        WRITE (*, *) nf_strerror(ierr)
663      END IF
664
665
666      ! write the 'times'
667      DO k = 1, 360
668#ifdef NC_DOUBLE
669        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
670#else
671        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
672#endif
673        IF (ierr/=nf_noerr) THEN
674          WRITE (*, *) 'writelim error with temps(k),k=', k
675          WRITE (*, *) nf_strerror(ierr)
676        END IF
677      END DO
678
679    END IF ! of if (is_master)
680
681    ! write the fields, after having collected them on master
682
683    CALL gather(phy_nat, phy_glo)
684    IF (is_master) THEN
685#ifdef NC_DOUBLE
686      ierr = nf_put_var_double(nid, id_nat, phy_glo)
687#else
688      ierr = nf_put_var_real(nid, id_nat, phy_glo)
689#endif
690      IF (ierr/=nf_noerr) THEN
691        WRITE (*, *) 'writelim error with phy_nat'
692        WRITE (*, *) nf_strerror(ierr)
693      END IF
694    END IF
695
696    CALL gather(phy_sst, phy_glo)
697    IF (is_master) THEN
698#ifdef NC_DOUBLE
699      ierr = nf_put_var_double(nid, id_sst, phy_glo)
700#else
701      ierr = nf_put_var_real(nid, id_sst, phy_glo)
702#endif
703      IF (ierr/=nf_noerr) THEN
704        WRITE (*, *) 'writelim error with phy_sst'
705        WRITE (*, *) nf_strerror(ierr)
706      END IF
707    END IF
708
709    CALL gather(phy_bil, phy_glo)
710    IF (is_master) THEN
711#ifdef NC_DOUBLE
712      ierr = nf_put_var_double(nid, id_bils, phy_glo)
713#else
714      ierr = nf_put_var_real(nid, id_bils, phy_glo)
715#endif
716      IF (ierr/=nf_noerr) THEN
717        WRITE (*, *) 'writelim error with phy_bil'
718        WRITE (*, *) nf_strerror(ierr)
719      END IF
720    END IF
721
722    CALL gather(phy_alb, phy_glo)
723    IF (is_master) THEN
724#ifdef NC_DOUBLE
725      ierr = nf_put_var_double(nid, id_alb, phy_glo)
726#else
727      ierr = nf_put_var_real(nid, id_alb, phy_glo)
728#endif
729      IF (ierr/=nf_noerr) THEN
730        WRITE (*, *) 'writelim error with phy_alb'
731        WRITE (*, *) nf_strerror(ierr)
732      END IF
733    END IF
734
735    CALL gather(phy_rug, phy_glo)
736    IF (is_master) THEN
737#ifdef NC_DOUBLE
738      ierr = nf_put_var_double(nid, id_rug, phy_glo)
739#else
740      ierr = nf_put_var_real(nid, id_rug, phy_glo)
741#endif
742      IF (ierr/=nf_noerr) THEN
743        WRITE (*, *) 'writelim error with phy_rug'
744        WRITE (*, *) nf_strerror(ierr)
745      END IF
746    END IF
747
748    CALL gather(phy_fter, phy_glo)
749    IF (is_master) THEN
750#ifdef NC_DOUBLE
751      ierr = nf_put_var_double(nid, id_fter, phy_glo)
752#else
753      ierr = nf_put_var_real(nid, id_fter, phy_glo)
754#endif
755      IF (ierr/=nf_noerr) THEN
756        WRITE (*, *) 'writelim error with phy_fter'
757        WRITE (*, *) nf_strerror(ierr)
758      END IF
759    END IF
760
761    CALL gather(phy_foce, phy_glo)
762    IF (is_master) THEN
763#ifdef NC_DOUBLE
764      ierr = nf_put_var_double(nid, id_foce, phy_glo)
765#else
766      ierr = nf_put_var_real(nid, id_foce, phy_glo)
767#endif
768      IF (ierr/=nf_noerr) THEN
769        WRITE (*, *) 'writelim error with phy_foce'
770        WRITE (*, *) nf_strerror(ierr)
771      END IF
772    END IF
773
774    CALL gather(phy_fsic, phy_glo)
775    IF (is_master) THEN
776#ifdef NC_DOUBLE
777      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
778#else
779      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
780#endif
781      IF (ierr/=nf_noerr) THEN
782        WRITE (*, *) 'writelim error with phy_fsic'
783        WRITE (*, *) nf_strerror(ierr)
784      END IF
785    END IF
786
787    CALL gather(phy_flic, phy_glo)
788    IF (is_master) THEN
789#ifdef NC_DOUBLE
790      ierr = nf_put_var_double(nid, id_flic, phy_glo)
791#else
792      ierr = nf_put_var_real(nid, id_flic, phy_glo)
793#endif
794      IF (ierr/=nf_noerr) THEN
795        WRITE (*, *) 'writelim error with phy_flic'
796        WRITE (*, *) nf_strerror(ierr)
797      END IF
798    END IF
799
800    ! close file:
801    IF (is_master) THEN
802      ierr = nf_close(nid)
803    END IF
804
805  END SUBROUTINE writelim
806
807  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
808
809  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
810    USE dimphy
811    IMPLICIT NONE
812
813    INTEGER nlon, type_profil, i, k, j
814    REAL :: rlatd(nlon), phy_sst(nlon, 360)
815    INTEGER imn, imx, amn, amx, kmn, kmx
816    INTEGER p, pplus, nlat_max
817    PARAMETER (nlat_max=72)
818    REAL x_anom_sst(nlat_max)
819    CHARACTER (LEN=20) :: modname='profil_sst'
820    CHARACTER (LEN=80) :: abort_message
821
822    IF (klon/=nlon) THEN
823       abort_message='probleme de dimensions dans profil_sst'
824       CALL abort_physic(modname,abort_message,1)
825    ENDIF
826    WRITE (*, *) ' profil_sst: type_profil=', type_profil
827    DO i = 1, 360
828      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
829
830      ! Rajout fbrlmd
831
832      IF (type_profil==1) THEN
833        ! Méthode 1 "Control" faible plateau à l'Equateur
834        DO j = 1, klon
835          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
836          ! PI/3=1.047197551
837          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
838            phy_sst(j, i) = 273.
839          END IF
840        END DO
841      END IF
842      IF (type_profil==2) THEN
843        ! Méthode 2 "Flat" fort plateau à l'Equateur
844        DO j = 1, klon
845          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
846          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
847            phy_sst(j, i) = 273.
848          END IF
849        END DO
850      END IF
851
852
853      IF (type_profil==3) THEN
854        ! Méthode 3 "Qobs" plateau réel à l'Equateur
855        DO j = 1, klon
856          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
857            rlatd(j))**4)
858          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
859            phy_sst(j, i) = 273.
860          END IF
861        END DO
862      END IF
863
864      IF (type_profil==4) THEN
865        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
866        DO j = 1, klon
867          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
868            rlatd(j))**4)
869          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
870            phy_sst(j, i) = 273.
871          END IF
872        END DO
873      END IF
874
875      IF (type_profil==5) THEN
876        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
877        DO j = 1, klon
878          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
879            *rlatd(j))**4)
880          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
881            phy_sst(j, i) = 273. + 2.
882          END IF
883
884        END DO
885      END IF
886
887      IF (type_profil==6) THEN
888        ! Méthode 6 "cst" valeur constante de SST
889        DO j = 1, klon
890          phy_sst(j, i) = 288.
891        END DO
892      END IF
893
894
895      IF (type_profil==7) THEN
896        ! Méthode 7 "cst" valeur constante de SST +2
897        DO j = 1, klon
898          phy_sst(j, i) = 288. + 2.
899        END DO
900      END IF
901
902      p = 0
903      IF (type_profil==8) THEN
904        ! Méthode 8 profil anomalies SST du modèle couplé AR4
905        DO j = 1, klon
906          IF (rlatd(j)==rlatd(j-1)) THEN
907            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
908              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
909            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
910              phy_sst(j, i) = 273. + x_anom_sst(pplus)
911            END IF
912          ELSE
913            p = p + 1
914            pplus = 73 - p
915            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
916              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
917            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
918              phy_sst(j, i) = 273. + x_anom_sst(pplus)
919            END IF
920            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
921          END IF
922        END DO
923      END IF
924
925      IF (type_profil==9) THEN
926        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
927        DO j = 1, klon
928          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
929            *rlatd(j))**4)
930          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
931            phy_sst(j, i) = 273. - 2.
932          END IF
933        END DO
934      END IF
935
936
937      IF (type_profil==10) THEN
938        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
939        DO j = 1, klon
940          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
941            *rlatd(j))**4)
942          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
943            phy_sst(j, i) = 273. + 4.
944          END IF
945        END DO
946      END IF
947
948      IF (type_profil==11) THEN
949        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
950        DO j = 1, klon
951          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
952            rlatd(j))**4)
953          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
954            phy_sst(j, i) = 273.
955          END IF
956        END DO
957      END IF
958
959      IF (type_profil==12) THEN
960        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
961        DO j = 1, klon
962          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
963            *rlatd(j))**4)
964          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
965            phy_sst(j, i) = 273. + 4.
966          END IF
967        END DO
968      END IF
969
970      IF (type_profil==13) THEN
971        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
972        DO j = 1, klon
973          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
974            rlatd(j))**4)
975          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
976            phy_sst(j, i) = 273.
977          END IF
978        END DO
979      END IF
980
981      IF (type_profil==14) THEN
982        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
983        DO j = 1, klon
984          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
985            *rlatd(j))**4)
986          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
987            phy_sst(j, i) = 273.
988          END IF
989        END DO
990      END IF
991
992      if (type_profil.EQ.20) then
993      print*,'Profile SST 20'
994!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
995
996      do j=1,klon
997        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
998      enddo
999      endif
1000
1001      if (type_profil.EQ.21) then
1002      print*,'Profile SST 21'
1003!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1004      do j=1,klon
1005        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
1006      enddo
1007      endif
1008
1009
1010
1011    END DO
1012
1013    ! IM beg : verif profil SST: phy_sst
1014    amn = min(phy_sst(1,1), 1000.)
1015    amx = max(phy_sst(1,1), -1000.)
1016    imn = 1
1017    kmn = 1
1018    imx = 1
1019    kmx = 1
1020    DO k = 1, 360
1021      DO i = 2, nlon
1022        IF (phy_sst(i,k)<amn) THEN
1023          amn = phy_sst(i, k)
1024          imn = i
1025          kmn = k
1026        END IF
1027        IF (phy_sst(i,k)>amx) THEN
1028          amx = phy_sst(i, k)
1029          imx = i
1030          kmx = k
1031        END IF
1032      END DO
1033    END DO
1034
1035    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
1036    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
1037    ! IM end : verif profil SST: phy_sst
1038
1039    RETURN
1040  END SUBROUTINE profil_sst
1041
1042END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.