source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyaqua_mod.F90 @ 3820

Last change on this file since 3820 was 3818, checked in by millour, 10 years ago

Some partial cleanup on uses of "dimensions.h" in physics.
At this point 3D gcm compiles and bench seems to run fine :-)
EM

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